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AN ENTROPY-BASED APPROACH TO NONLINEAR STABILITY 


Marshal L. Merriam 
ABSTRACT 


Many numerical methods used in Computational Fluid Dynamics (CFD) incorpo- 
rate an artificial dissipation term to supress spurious oscillations and control nonlinear 
instabilities. The same effect can be accomplished by using upwind techniques, sometimes 
augmented with limiters to form Total Variation Diminishing (TVD) schemes. An analysis 
based on numerical satisfaction of the second law of thermodynamics allows many such 
methods to be compared and improved upon. For example, certain TVD schemes tend 
to “square” a smooth pulse. These can be detected a priori by their negative entropy 
production rates. 

A nonlinear stability proof is given for discrete scalar equations arising from a conser- 
vation law. Solutions to such equations are bounded in the £ 2 norm if the second law of 
thermodynamics is satisfied in a global sense over a periodic domain. It is conjectured that 
an analogous statement is true for discrete equations arising from systems of conservation 
laws. 


Stability in the L 2 norm is not sufficient to exclude expansion shocks, oscillations, 
and other unphysical phenomena. Numerical experiments suggest that a more restrictive 
condition, a positive entropy production rate in each cell, is both necessary and sufficient 
to exclude such phenomena. 

Construction of schemes which satisfy this condition is demonstrated for linear and 
nonlinear wave equations and for the one-dimensional Euler equations. For the linear wave 
equation all schemes which satisfy a cell entropy inequality are formally TVD. The converse 
is not true. The scheme for the Euler equations makes use of a remarkable identity to avoid 
the usual frozen coefficient approximation. Results are shown. 

Since this form of dissipation is based on classical thermodynamics, it extends natu- 
rally to all types of fluid dynamics problems. In particular, artificial dissipation require- 
ments for the Navier-Stokes equations, compressible and incompressible, may be quanti- 
tatively assessed. 

Most existing forms of dissipation (including those mentioned above), treat problems 
in two and three dimensions as a sequence of one-dimensional operators. The present work 
does not require this approximation, and therefore extends more naturally to unstructured 
meshes. It holds promise for use with multigrid and grid refinement methods. 


Chapter 1. INTRODUCTION 


1.1 The Importance of Nonlinear Stability 

For some time now it has been possible to obtain numerical solutions to physical prob- 
lems described by systems of conservation laws. Among these problems is heat conduction, 
e.g., the determination of temperature distribution in a solid heated at a boundary. Exist- 
ing numerical solution techniques for this problem are satisfactory in every respect. That 
is: they are feist, accurate, and robust. The same cannot be said of numerical solution 
techniques for problems in compressible fluid mechanics that are governed by the Navier- 
Stokes equations. While good comparison with experiment can sometimes be obtained, 
good error estimates usually cannot be. 

Recent advances in Computational Fluid Dynamics (CFD) have made it possible to 
obtain reasonable solutions to the Navier-Stokes equations. Such solutions are generally 
sensitive to the details of the computational grid, to the time step, to the turbulence 
model, to the boundary conditions, and to the details of the artificial dissipation used. All 
of these are adjusted in practice to improve the agreement between the numerical solution 
and a corresponding experimental measurement or analytical result. The large number of 
adjustments and the lack of a good error measure make it hard to use CFD in cases where 
there are no experiments. 

Navier-Stokes solvers tend to lack robustness. The more complicated the solver, the 
trickier it becomes to keep the computation stable and accurate. The usual scapegoats, 
cited with about equal frequency, are the grid, the time step, the turbulence model, the 
boundary conditions, and the artificial dissipation. In other words, selection of the proper 
adjustable parameters is critical. In the absence of theoretical guidance, a trial-and-error 
procedure is generally required to get a good answer. 

Finally, there are questions of generality. There is considerable current interest in 
Navier-Stokes solutions involving chemically reacting flows about complex geometries. Al- 
gorithms which rely on a characteristic diagonalization are unsuitable when the character- 
istics are not known, real, or computable, as will surely be the case for chemically reacting 
flows. Similarly, algorithms which rely on dimensional splitting are difficult to apply to 
the unstructured grids being used to model complex geometries. 

Numerical methods require speed, accuracy, and robustness to be useful. Faster com- 
puters improve the effective speed of a given method. Larger memories allow more grid 
points, thereby improving accuracy. Robustness, on the other hand, cannot be achieved 
through improvements in computer hardware. The aim of this work is to address the 
question of robustness. The approach explored here provides a physical, quantitative ex- 
planation for artificial viscosity. The improved understanding this requires provides insight 
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to other aspects of CFD. For example, the effects of time advance on stability are explored 
in Chapter 4 and a way of judging grid quality is introduced in Chapter 10. 

1.2 Early Work in Computational Fluid Dynamics 

The first attempts at numerical simulation of nonlinear convection were not completely 
successful. The solutions tended to develop numerical oscillations which grew exponen- 
tially as the solutions progressed in time. Approximating spatial derivatives proved to be 
especially troublesome; the most accurate formulae proved to be the least stable. What 
was lacking at the time was a physical and mathematical understanding of the convection 
process. 

Much of this confusion was alleviated by an early monograph (ref. 1) in which Peter 
Lax detailed the theory of shock waves and sonic expansions. Among the ideas presented 
were the need for an entropy condition to make the solution unique, and the use of artificial 
viscosity to enforce an entropy condition. Domain of dependence ideas were introduced 
within the framework of characteristic theory. At the end of the work are some simple 
numerical methods and proofs that they lead to an accurate answer (in the limit of very 
fine grids). The ideas presented in Lax’s monograph have been very important in the 
development of CFD. Their influence is apparent in every major scheme in use today, and 
in the work presented here. 

An early success was the work of Murman and Cole (ref. 2) on the transonic small 
disturbance equations. The main advance in this work wits the Murman-Cole switch which 
used domain of dependence ideas, cell-by-cell, to choose between two approximations for 
the streamwise derivatives. One stencil would be used in subsonic flow, another in su- 
personic flow. Though limited to transonic or subsonic flow about slender bodies with 
attached flow, this work represented a tremendous advance over the calculations that pre- 
ceded it. To quote Robert MacCormack (ref. 3) “ The use of these methods spread rapidly 
and airfoil sections were no longer designed experimentally by cutting and filing.” 

In the mid-seventies, the first solutions to the Navier-Stokes equations appeared. 
Central to this advance was the work of Steger (ref. 4) and Pulliam and Steger (ref. 5) 
in their ARC2D and ARC3D codes. In addition to the use of generalized coordinates 
(but not generalized topologies), these codes used implicit time advance and an artificial 
dissipation based on fourth differences. They also used the Baldwin-Lomax turbulence 
model (ref. 6), a simple, empirical model that replaced turbulence with an equivalent 
eddy viscosity. Although suffering from oscillations near shocks and having very poor 
convergence by today’s standards, these codes represented a revolutionary advance in the 
ability to calculate separated flows and model high Reynolds number flows. 

In the late seventies and early eighties, algorithms began to appear (refs. 7,8, and 9) 
which incorporated a characteristic decomposition of the inviscid fluxes. It was noted 
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that the flux Jacobian with respect to the conservative variables had real eigenvalues and 
eigenvectors which were not too difficult to compute. This allowed the use of domain of 
dependence ideas in the Euler equations. The algorithms which resulted had no need for 
artificial dissipation of the type used in ARC2D as they were naturally dissipative. This 
feature removed the need to choose the coefficient for artificial dissipation. It also removed 
any control over it. In spite of the major increase in physical understanding included 
in these algorithms, they continued to exhibit minor numerical oscillations which had no 
physical basis. 

To get rid of these annoying oscillations, researchers attacked them head on, produc- 
ing schemes which, at least for scalar equations in one dimension, were guaranteed not to 
have any wiggles. The first scheme of this type, Flux Corrected Transport (ref. 10) demon- 
strated the concept. Later the concept of Total Variation Diminishing (TVD) schemes was 
formalized by Harten (ref. 11) and implemented for gasdynamics. Although Lax referred 
to proofs (refs. 12 and 13) that the desired physical solution has the TVD property, it 
does not follow that all TVD schemes converge to a physical solution. Indeed, schemes 
have been introduced (ref. 2) which satisfy the TVD condition and yet support unphysical 
expansion shocks. In the rush to kill off those annoying wiggles, some researchers have 
ignored the entropy inequality that is fundamental to getting a correct solution. 

In this decade, some research has been done on schemes which satisfy an entropy 
inequality. Some of this work, done by Osher (ref. 14) and Tadmor (ref. 15), will be 
explored in later chapters. It concerns a numerical approximation to the entropy inequality 
on a cell-by-cell basis. From a different perspective, Hughes, Franca, and Mallet (ref. 16) 
have concentrated on a finite element method which guarantees global satisfaction of the 
second law of thermodynamics in the steady state. The local oscillations which still occur 
offer evidence of local violations of the second law, even though it is globally satisfied. The 
thrust of the work presented here will be to satisfy the second law of thermodynamics on 
a cell-by-cell basis. 

1.3 Thesis Organization and Summary of Results 

The goal of this research is to analyze the requirements for nonlinear stability and 
produce an algorithm to show that these can be met. Chapter 2 establishes notation 
and shows some important properties of partial differential equations (PDE’s) which come 
from conservation laws. These properties are exploited to show that solutions of scalar 
hyperbolic equations in multidimensions are bounded for all time. For systems of equations, 
analysis shows only that the total entropy in the domain is bounded; the question of bounds 
on the solution remains open. 

In Chapter 3 the nonlinear stability property of the PDE is extended to the semi- 
discrete difference equations. The semi-discrete analog of the second law is derived. Chap- 
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ter 3 goes on to discuss, in general terms, how to construct space differencing schemes 
which satisfy a cell entropy inequality. Such schemes can be constructed with at least 
second order accuracy. 

In Chapter 4 the nonlinear stability property of the PDE is extended to the fully dis- 
crete difference equations. It is shown that satisfying a semi-discrete cell entropy inequality 
is sufficient to satisfy a fully discrete one; provided that implicit Euler time advance is used. 
A modified Crank-Nicolson time advance for scalar equations also has this property. 

In Chapters 5 through 7, examples are given for problems that increase in complexity 
up to the quasi-one-dimensional Euler equations. These examples show how the ideas 
given in Chapters 2, 3, and 4 can be applied in practice. Chapters 8 and 9 show how these 
new ideas relate to existing schemes. 

In Chapter 10 a good deal of future work is outlined. A scheme is not very useful 
if it has limited applicability. On the other hand, it is not worthwhile to extend an 
impractical scheme; better to clean up the simple case first. For this reason, Chapter 10 
is limited to describing in some detail how these extensions might be made and discussing 
the implications of these ideas on several of the major problems that still confront us in 
CFD. The author feels that these ideas constitute a powerful new tool, making possible a 
whole new approach to these problems. Some obvious approaches will be outlined. Finally, 
Chapter 11 provides a summitry. 
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Chapter 2. STABILITY AND THE SECOND LAW 


2.1 A Physical Description of Entropy 

The Euler equations for gasdynamics in one dimension have a well-defined conser- 
vation form and do not allow heat transfer. They also have no viscous terms. In one 
dimension the Euler equations can be written 


9>t + f>z = 0 


( 2 . 1 . 1 ) 


where 



’ P " 


pu 

9 = 

pu 

and / = 

pu 2 +p 


e 


.(e + p)u 


( 2 . 1 . 2 ) 


The quantities p and u represent the local mass density and the local velocity respectively. 
The quantity e represents the total energy density. The total energy can be expressed as 
the sum of the thermal and kinetic energies of the gas. 


e = pc v T + -pu 1 


(2.1.3) 


The variable T represents the local temperature and c v represents the heat capacity at 
constant volume. For a thermally and calorically perfect gas, p — pRT. This equation of 
state, and the relations 7 = ^ and R = cp — c», allow derivation of an expression for 
pressure. 


P = (7 - l)(e - j/”* 2 ) 


(2.1.4) 


Equations (2.1.1) are formally correct only as long as the derivatives are defined. The 
corresponding integral equations are more lenient. These are 


f q[t + At] dv - f q[t\ dv + f <f 
Jv Jv Jt Ja 


f • n dA dr = 0 


av 


(2.1.5) 


where n represents the outward facing unit normal to the surface dV (in one dimension 
the surface integral reduces to a difference). To derive (2.1.1) from (2.1.5), the divergence 
theorem is used to convert the surface integrals into volume integrals. The limit is then 
taken as the space time control volume becomes infinitesimally small. At this limit, the 
integrand must vanish if the integral vanishes. Such a step is valid only in the absence 
of a discontinuity. Because of this caveat, the use of (2.1.1) is prone to difficulties when 
applied to flows with discontinuities. 


Equation (2.1.5) represents three separate conservation laws, one for each component 
of q. The first equation corresponds to conservation of mass, the second corresponds 
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to conservation of momentum, and the third equation is a statement of the first law of 
thermodynamics. 


« 





(e + p)u • n dA dr = 0 


( 2 . 1 . 6 ) 


The conserved quantity in (2.1.6) is the total energy, the sum of the internal energy 
density ( pc v T ) and the kinetic energy density (|/m 2 ) integrated over any control volume. 
These two energy forms are not (and should not be) individually conserved; the first law 
allows arbitrary transfers of energy between them. 

A fluid which is solely governed by the conservation laws in (2.1.5) can display very 
unusual behavior. For example, it can transfer all of its internal energy to kinetic energy. 
The result is a very fast, very cold jet of air. This Carnot air conditioner is allowed under 
conservation of mass, momentum and energy, yet is never observed in practice. To explain 
this discrepancy requires introducing the second law of thermodynamics. Under the second 
law, certain transfers between interned energy and kinetic energy are not possible. 

Three independent variables are required to describe the fluid of (2.1.5). A funda- 
mental result of thermodynamics suggests that such a gas has two reversible work modes 
(ref. 17). These are known for the one-dimensional Euler equations. If work is expended 
to increase fluid kinetic energy at constant pressure and temperature, it can be completely 
recovered by decelerating the flow. This is called an isentropic acceleration and decelera- 
tion. If work is expended to compress the fluid isentropically at constant velocity (raising 
its temperature in the process), it can be completely recovered in a subsequent expansion 
(which lowers the temperature). Through these processes, the kinetic energy and internal 
energy (temperature) can be independently and reversibly varied. Other reversible state 
changes are a combination of these two processes. For example, the work required for an 
isentropic compression can be provided by an isentropic deceleration. Similarly, the work 
required for an isentropic acceleration can be provided by an isentropic expansion. 

In addition to reversible state changes, there are irreversible ones. For example, 
heat transfer across a finite temperature difference is irreversible as is viscous momentum 
transfer. Although not explicitly included in (2.1.5), these mechanisms are significant 
within the profile of shocks, the main source of irreversibility in the Euler equations. 
The word irreversible implies that certain state changes (like transferring heat against 
a temperature gradient) are not possible; they are never observed in the laboratory on 
macroscopic scales. The Carnot air conditioner described above requires the use of state 
changes that are impossible in this sense. 
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The second law of thermodynamics provides a way of discriminating between the 
possible state changes and the impossible. It can be written in a form similar to (2.1.5) 



dv dr = 


[ S[t + At] dv - f 5[t] dv + f (f F • n dA dr > 0 
Jv Jv Jt Jev 


(2.1.7) 


In this equation, P s represents the rate of production of entropy per unit volume, 5 is the 
entropy per unit volume, and F is the entropy flux per unit area. For reversible processes, 
p a = 0. For irreversible processes P t > 0. Processes for which P t < 0 are never observed. 


The quantities S and F depend on the fluid being studied and on the idealizations 
being applied. For the Euler equations in the absence of heat transfer, these are 


S = ps 
F = pus 

where s is the nondimensional specific entropy. 


( 2 . 1 . 8 ) 

(2.1.9) 


From Reynolds (ref. 18), comes the expression for the entropy difference between an 
arbitrary state, q' and a reference state g< >, for a perfect gas. 


- 50 


c v 



( 2 . 1 . 10 ) 


Equation (2.1.10) suggests nondimensional forms for entropy, density and pressure. Using 
these gives 

»=log(£) (2.1.11) 

which is used throughout this work. The second law (2.1.7) is not limited to the Euler 
equations or to an ideal gas. It applies to all of continuum mechanics. 

2.2 A Mathematical Description of Entropy 


In addition to the physical view of entropy given above, there is a well-developed 
mathematical theory of entropy (from Lax (ref. 19) and Krushkov (ref. 20)) given in the 
monograph previously mentioned (ref. 1). This theory, now well known, holds that any 
entropy pair ( S , F ) must have two properties. 

S, qq < 0 ( 2 . 2 . 1 ) 

F„ = S, q f„ ( 2 . 2 . 2 ) 

where S is the the entropy per unit volume and F is the entropy flux. 
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The quantities which appear in (2.2.1) and (2.2.2) tire defined as follows 


S 


m 


S 


>« 


’ d 2 s 

d 2 s 

d 2 s * 


dfi 

dh 

dh] 

dqidqi 

dqidq 2 

dq\dq n 


dqi 

dq 2 

dq n 

d 2 s 

d 2 s 

d 2 s 


dh 

dh 

Oh 

dqidqi 

dq 2 dq 2 

dq 2 dq n 

and /„ = 

dq\ 

dq 2 

d g„ 

• 

d 2 S 

d 2 s 

d 2 s 


dfn 

dfn 

dfn 

. dq n dqi 

dq n dq 2 

&Qn &Qn m 


. dqi 

dq 2 

dq n . 


' ds 

dS 

as 

• • • 

and jpjg ~ 

dF 

dF 

dF 

• • • > ■ 

.dqi' 

dq 2 * 

’ dq n \ 


.dqi' 

dq 2 

’ dqn J 


(2.2.3) 


(2.2.4) 


The convexity condition (2.2.1) forces irreversible processes to run in the correct direc- 
tion and produce entropy. A consequence of the compatability condition (2.2.2) is that 
reversible processes do not produce entropy. The following two examples illustrate these 
effects for the Euler equations. 


It is well known that if a quantity of fluid is placed in an isolated box, the equilibrium 
condition will be the one for which entropy is a maximum. This occurs when till of the 
fluid is at the same state. This state, g, is simply the constant state with the correct total 
mass, momentum and energy for the box. 


9 = 



(2.2.5) 


It is the nature of isolated boxes that there is no flux through their boundaries. This, 
combined with the conservation law (2.1.5), is sufficient to guarantee that q doesn’t change 
with time. The state transition from g(x) to q is physically irreversible so the entropy 
contained within the volume should increase. This property follows mathematically if the 
convexity condition is met, as can be shown by using Taylor’s theorem with remainder as 
follows. 

S = S + S„(g - q) + i(q - ?) T S,**(g - q ) (2.2.6) 

Here the quantity q is used to define 5 = 5(g) and S, q = 5, ? (g). The term 5,** is evaluated 
at some unspecified state, g**, which makes (2.2.6) an equality. The existence of such a 
state is guaranteed by the intermediate value theorem. 

Both sides of (2.2.6) can be integrated over a control volume (which encloses the box) 
to give 

Jv Sdv = Jv S dv + S* q jy ( ? ” *) dv + \ J v ( ? “ - q) dv ( 2 . 2 . 7 ) 
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The second term on the right-hand side vanishes by the definition of q. 

h SdV = iv SiV *\Jv (q - 5)TS '« (? ~ ’ )<f ” ( 2 - 2 - 8 ) 

The integrand of the last term is a proper quadratic form since S,** is a negative , 

definite matrix. This implies that the last term is nonpositive. In this way, convexity is 
used to show that the maximum entropy occurs when the equality holds, i.e., when q = q 
over the entire cell. < 

The compatability condition (2.2.2) addresses reversible processes. The Elder equa- 
tions (2.1.5) do not contain terms for heat transfer, viscosity or any other dissipation 
mechanism. As long as the solution is smooth, the flow processes are reversible and no 
entropy production should occur. Such a statement follows mathematically if the com- 
patability condition is met. This is shown by first multiplying (2.1.1) on the left by 5,,. 

"H Syqffx — 0 (2.2.9) 

In smooth regions (where the derivatives are finite) the chain rule can be applied as follows 

S,t + = 0 (2.2.10) 

At this point we assume that (2.2.2) holds and make the substitution. 

S, t +F, q q, t = 0 (2.2.11) 

Using the chain rule one more time leads to 

S,t + F,* = 0 (2.2.12) 

which is the differential statement of conservation of entropy. This shows that if (2.2.2) 
holds and if the solution is smooth, then entropy will be conserved and P t = 0. Such a 
statement is known to be true for the Euler equations (ref. 21), in which entropy production 
is nonzero only at discontinuities. 

The definitions of S and F are not always unique. For example, Harten (ref. 22) 
showed a family of entropy pairs for the Euler equations. Furthermore (again from Lax 
(ref. 1)) if q contains more than two independent variables, (2.2.2) represents a system of * 

partial differential equations that is overdetermined and may have no solution. Fortunately, 
for the specific case of the Navier-Stokes equations in three dimensions, Mallet (ref. 23) 
showed that S and F exist and are essentially unique. 
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2.3 A Relationship Between Mathematical and Physical Entropy 


At this point, a question arises: Does the mathematical entropy described in Section 
2.2 correspond to the physical entropy described in Section 2.1? At least for some impor- 
tant cases, the answer is yes. For the one-dimensional Euler equations, q and / are defined 
as in (2.1.2) 



’ P 


pu 

9 = 

pu 

/ = 

pu 1 +p 


e 


(e + p) u 


(2.3.1) 


As mentioned in Section 2.1, the entropy pair (S,F) suggested by thermodynamics is 


where 


S = ps 

(2.3.2) 

F = pus 

(2.3.3) 


(2.3.4) 


Given these definitions for q , /, 5, and F, it is possible to check the physical entropy 
and entropy flux against the mathematical constraints (2.2.1) and (2.2.2). The first check 
is whether S m < 0. The derivatives are easily carried out using the symbolic algebra 
program MACSYMA; the matrix is 



4 


7 P. , . 
(7-») a P 


2 

pu 2 — e 


*£ + e -pu 
-pu p 


The determinant of S, qq is also easily computed by MACSYMA and is 




(7 ~ l) 4 
P 3 


(2.3.5) 


(2.3.6) 


which indicates that S m is nonsingular, at least for p > 0. Furthermore there exists, for 
the Euler equations, a matrix Y' 1 such that 


S m = ~(Y') t (Y-') 


(2.3.7) 


This implies that S, qq is negative semi-definite. Since (2.3.6) rules out the singular case, 
S tqq is, in fact, negative definite, demonstrating that (2.2.1) is satisfied when S is given by 
(2.3.2). 
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The matrices Y and Y' 1 are given below. For convenience, the speed of sound c = 
has been introduced. 


Y = 


Pi P 2 Pi 

Piu p 3 (u - c) p 2 (u + c) 

*(£+•“+ r£r) 


(2.3.8) 


-=(*) 


r ^ c 3 ~(7~^)^ 2 2(7— l)tt —2(7—1) 

Pi Pi Pi 

(7— l)tt 2 , uc -(y-l)u-c (7I) 

2p2 p2 P 2 p2 

(7-1)1** _ tic -(7-l)«+c (t- 1) 

2 p2 p2 p2 P2 


A = V /f’ ft = / 


27(7 - 1) 


(2.3.9) 


(2.3.10) 


Identity (2.3.7) has been verified by construction using MACSYMA. This identity is a new 
result. It is discussed further in Chapter 7. 

The second requirement on the mathematical entropy is the existence of an entropy 
flux F such that 

F„ = S„f„ (2.3.11) 

Such an entropy flux is given by (2.3.3). This is easily verified with MACSYMA. For 
reference the individual terms are given here. 


F w = [i2^e5!-7u, telle!, teii e ] 
S., = [ < -(7 + l) + ^, -lac^ies, fa=Ue ] 


/»« “ 




1 

-(7 “ 3)« 


L( 7 -1)u 3 -^ 5-|( 7 -I )* 1 


0 

(7-1) 

ju 


(2.3.12) 

(2.3.13) 

(2.3.14) 


Since satisfaction of both requirements has been verified, (2.3.2) and (2.3.3) define a 
valid entropy pair for the Euler equations in the mathematical sense. It is well known that 
this entropy pair also describes the physical entropy density and entropy flux (ref. 23). 

This example, and those of the previous section provide indications that the math- 
ematical and physical entropy statements are compatible; if one is true then the other is 
true. Such a conclusion is not new; these examples are included for the benefit of readers 
unfamiliar with this somewhat esoteric subject. 


1 
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2.4 The Sign Convention for Entropy 


When entropy was first explored from a physical point of view by Carnot and others, 
the entropy statement was viewed almost theologically. The entropy of the universe is 
increasing; everything goes toward chaos. Given this point of view, it seems natural that 
entropy production should be a positive quantity. For scalar conservation laws with a 
single dependent variable u, S = — u 2 . The second law requires that this increase with 
time for any closed system. 

By contrast, the first mathematical treatments of conservation laws were concerned 
with bounds on the solution. The natural statement of such bounds for scalar conservation 
laws is that u 2 is conserved or decreases with time. This seemed more natural than an 
equivalent statement about -u 2 . When the concept of entropy was later introduced, con- 
sistency of notation demanded that entropy should decrease with time. As a consequence, 
many papers on the subject consider mathematical entropy to be a decreasing quantity, 
positive entropy production rates to be impossible and S, qq to be positive definite. 

This difference in sign convention is pointed out to help the interested reader avoid 
confusion when reviewing the literature on this subject. There is no substantive difference 
between the two points of view. As a purely arbitrary choice, this work adheres to the 
physical convention. 


2.5 Nonlinear Stability 

Nonlinear stability is the holy grail of numerical analysis. Always sought, occasion- 
ally sighted, this goal has eluded researchers for many years. In this study, a promising 
approach is examined in some detail. Some of the ideas in this approach are due to Dutt 
(ref. 24) who several years ago offered a stability proof subject to certain questionable 
assumptions which will be pointed out along the way. What follows is a greatly simplified 
and shortened version of his proof, following the same general outline. Though a complete 
proof still eludes us, the goal seems much closer than it was. 

The basic idea is to first show that the total entropy in the domain is bounded. This 
can be done by showing that it steadily increases and that it cannot increase above a 
certain level. The next step is to show that if the entropy is bounded, the solution is 
bounded in a certain norm. 

Conjecture. 

Over some domain which encloses volume fl and is enclosed by a surface dfl 
Given: 

<p f • n dA = 0 

Ja n 


(2.5.1) 
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/ «[*) rfv - 

Ja 

f q[0]dV+ ( / f • n dA dr = 0 
J o Jo J an 

(2.5.2) 



s'®'' 

• 

s 

£ 

ii 

o 

(2.5.3) 


f S[t]dV- 
Ja 

/ 5[0] dV + / <f F ‘XidAdr >Q 
Ja Jo Jaa 

(2.5.4) 

and defining 


1 t 




m=fij a mdv 

(2.5.5) 

There exists a 

value q * which can 

be determined from the initial conditions q[0] such that 


o<-i / (t-t) T sr„(,-i)dv< sm)-sm)dv< I sm)-sm)w 

2 J a J n J n 

(2.5.6) 


Notice that -± / n (q - q) T S>* q (q - q) dV is a norm, since 5,*, is negative definite. It is in 
this norm that the solution remains bounded. 

REASONING: 

Combining the conservation law (2.5.2) with the periodicity constraint (2.5.1) gives 
the relation 

f q[t)dV= f q[0]dV (2.5.7) 

J n Ja 

which says that the total amount of q in the domain is a constant over time. Dividing by 
the domain volume and applying definition (2.5.5) to both sides gives 

?[i] = »[0] (2.5.8) 

which is to say that q is also a constant over time. This being the case, the argument can 
be dropped and <j[f] can be expressed simply as q. 

At this point an earlier result (2.2.8) can be invoked using the entire domain as the 
control volume. 

/ S[t] dV= f SdV + l f (q- q) T S”(q - q) dV (2.5.9) 

J n Ja 1 Ja 

As before S = S(q). Also, the quadratic term is negative at each point in the domain so 


(2.5.11) 


Combining (2.5.4) with (2.5.3) implies that 


f S[t]dV> f S[0]dV 
J n Jo 


which implies that the total amount of entropy in the domain is a nondecreasing function 
of time. Combining (2.5.11) with (2.5.10) gives 


f SdV > f S[t]dV> f S[0]dF 
Jo Jo Jo 


(2.5.12) 


so the total entropy in the domain is bounded for all time. Note that both bounds can be 
computed from the initial conditions. Subtracting the upper bound from each term and 
multiplying by —1 gives 


f (S-S[0])dV> f (S-S[t))dV> 0 
Jo Jo 


(2.5.13) 


In Dutt’s proof, this is an intermediate result. There is a significant difference, how- 
ever, in that he uses a particular value of q which does not depend solely on the initial 
conditions. Rather, he limits his proof to the Navier-Stokes equations and assumes bounds 
on the density and temperature for all time. This results in a less precise bound and, more 
importantly, assumes a good deal of what he is trying to prove. 

This concludes the first portion of the proof. Equation (2.5.13) guarantees a bound 
on / n (S — 5[t]) dV. The second portion of the proof deals with the stability implications 
of this statement. 

Since the first portion of the proof dealt with domain integrated quantities, a suitable 
stability definition is 

IkWII < *ll«(0)ll (2.5.14) 

for 0 < t < T. Any norm can be used, including the Li norm. It is assumed that the 
scheme integrates up to a fixed time T > <o- It would be nice if k was near unity instead 
of some large value, but this is not required. 

The middle term of (2.5.13) can be calculated using (2.5.9). 

f (S^SlO))dV> [ (S-S\t))iV = -l f (q-q) T S"(q-q)dV >0 (2.5.15) 

J n J n n 


Notice that — i / n (q — q) T S,**(q — q)dV is almost a suitable norm. It is always pos- 
itive, since the integrand is a proper quadratic form and 5,** is negative definite. The 
problem is that q** is not determined. Furthermore it is different for each control volume. 
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(2.5.16) 


These problems can be overcome if a single q * exists such that 

-1)W>-\j a (i- «) T s>;,(« -q)dv> o 

where S, qq = S, qq (q*). The state q* must not be a function of time. Therefore it must 
depend only on the initial conditions. If such a state can be found, then a suitable norm is 

II# = J a (q -q)dV (2.5.17) 

Such a norm is bounded for all time by the relation 

f ($-S[0])iV> f (S-S(())dV = -i f (q - q) T S£(q ~q)dV > ||g|| 3 > 0 
J n J n * J n 

(2.5.18) 

It is not obvious how to find q * in general, which is why this is a conjecture and 
not a theorem. Some comments are in order, however. First, (2.2.6) is exact. A similar 
statement which is only approximately true is 

S(q) «3-5„(9-?) + i(i- q) T S m (q - ?) (2.5.19) 

This comes from using a truncated Taylor series instead of Taylor’s theorem with remainder 
so S,** « S m , with the approximation improving steadily as q[t] approaches q. This 
suggests that q might meet the requirements for q*. It can be determined from the initial 
conditions and 

-\J(q-lfS,„(q-q)dV>0 (2.5.20) 

The remaining requirement, (2.5.16), may be true, but is not easily proved in general. 

A second comment is that — S',** is positive definite. The evaluation state q** varies 
over the time and space domains of the problem. Presumably there is then some minimum 
value of S',*’. (This is tacitly assumed by Dutt in his proof.) The value of q at which this 
occurs could then be used for q*. Such a value of q* meets the requirements of (2.5.16) 
but is not a function of the initial conditions. Such an a posteriori norm has no predictive 
value. This is because it may assign a very small weight to a mode which is growing 
exponentially while assigning larger weights to modes which decay. While such a norm 
may decay for a time, it will eventually grow without bound. 

In summary then, this work differs from Dutt’s in everything but the general outline. 
It removes some unnecessary assumptions and an unjustified one. It replaces general 
bounds with specific ones. Most important, the analysis is more generally applicable and 
easier to follow than the pioneering work of Dutt. An unfortunate effect of all this progress 
is the downgrading of the result from a proof to a conjecture, a consequence of removing 
the assumption that q * exists in general. 
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Scalar Equations 

There is considerable work in the literature concerning scalar PDEs of the form 


ti,t + /,» = 0 

(2.5.21) 

where u is traditionally used in place of q in this context. It is easy to construct a valid 
entropy pair for such equations. The entropy 

1 

III 

<0 

(2.5.22) 

This certainly meets the convexity requirement since 


S ,vti = 2 0 

(2.5.23) 

Finally the compatability constraint 


F) U = Syvfiv 

(2.5.24) 

can be used to define F , if f(u) is differentiable. The quantities S, u 
constructed, multiplied, and integrated to give F(u). 

and /, u can be 

The importance of this rather obvious result is that the above conjecture becomes a 
theorem for scalar equations. This follows because S iUU is a constant, independent of the 
evaluation point. The result is that for scalars 

/($-S[0])<fl'>||( ? — 5)111 >0 

Jo 

(2.5.25) 

In fact, using (2.5.22) and (2.5.12) gives 


Ml > llff(*)lls > umi 

(2.5.26) 


Using the L 2 norm and k = 1 , (2.5.14) holds. Thus, the second law is sufficient for £2 
stability of all scalar hyperbolic PDEs. This result holds for multidimensional periodic 
domains and is probably extendable to unbounded domains as well. 


Euler Equations 

In the case of the Euler equations in one dimension, the conjecture has an interesting 
interpretation. Recall from Section 2.3 that for these equations 

Sm = -rffr 1 ) (2.5.27) 

If the conjecture is true, q * exists and one can construct new variables of the form 

q = Y-'q (2.5.28) 
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where the matrix Y' 1 is always evaluated at q*. The average values, q, can also be subjected 
to this transformation. If q* exists the result can be written 

jf (S - S[0])<(V > l||(,~?)III > 0 (2.5.29) 

Thus, the second law is sufficient for £2 stability of the one-dimensional Euler equations 
in the transformed variables, provided that the conjecture holds. 

The author has made some effort to establish the existence of q* in this case. In 
particular, an attempt was made to show that 

S,’„ - S " > 0 (2.5.30) 

which would guarantee the inequality (2.5.16). This required bounds on q* * and ultimately, 
bounds on q. Thus the attempt failed due to circular reasoning. The conjecture that q 
is bounded holds only if q can be shown to be bounded by some independent means. It 
seems likely that (2.5.16) will have to be shown in an integral sense and not in a pointwise 
sense if this conjecture is to be upgraded to a proof. 

In the remainder of this work the conjecture (2.5.6) is assutned to hold. This means 
that the satisfaction of a global entropy inequality implies nonlinear stability. The satisfac- 
tion of a cell-by-cell entropy inequality is certainly sufficient to satisfy a global inequality 
as well. 
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Chapter 3. THE SEMI-DISCRETE CASE 


The previous chapter discussed conditions for L 2 bounds on the solution of systems 
of hyperbolic conservation laws. These conditions and the equations themselves assume 
an arbitrary distribution of q over the domain. Such a distribution requires an infinite 
amount of information to express, a feature that is not conducive to computational solution 
techniques. To avoid such a problem, this chapter will derive equations which are discrete 
in space but continuous in time. These are called the semi-discrete equations. Later, it will 
discuss a semi-discrete form of the second law and show that the conjecture of the previous 
chapter is not affected by the semi-discrete approximation. Finally, it will explore ways 
to construct numerical schemes which satisfy the semi-discrete equations while meeting 
semi-discrete constraints. 


3.1 Deriving the Semi-Discrete Equations 

As a first step, assume that the domain of interest is tessellated into a finite number 
of control volumes, each of which have a finite size. These can be identified by a single 
index j. Furthermore the surface bounding the j th volume is denoted dVj. 


The integral equation describing such a volume is 

fi+At f 


f q[t + At]dv — f q[t]dv+ f <£ f - ndAdr = 0 (3.1.1) 

JVj JVj Jt JaVj 


Because each control volume has a finite extent, no singularities can occur in the limit as 
At — ► 0. Proceeding to this limit and dividing by At, 


d_ 

dt 



f • n dA = 0 


(3.1.2) 


The time indexes have been dropped because (3.1.2) holds at each time. At this point, the 
equations can be simplified by a definition. 


<U = 



(3.1.3) 


In one dimension there is a particularly simple numbering in which adjacent cells have 
adjacent numbers as shown in figure 3.1. 


j* Control Volume 



X j-1 X j-'A X j X j*'A X j + ] 

Figure 3.1. — A one-dimensional domain. 
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The surface integrals in one dimension are easily computed. Let the surface separating 
the j th cell and the j + l <fc cell be denoted j -j- so that dVj is composed of the two points 
Xj + i and Xj_i. Finally, let / = f • i, where i is the unit vector in the * direction. This 
allows the surface integrals to be carried out in closed form. For example 

jf f-n dA = f j+ ^-f M (3.1.4) 

In summary, the integral equations in one dimension reduce to 

v j(9j)>t + fj+$ ~ fj-\ ~ ® (3.1.5) 

Here the comma notation has been used to indicate differentiation. Traditionally the 
volume is divided out. In one dimension it seems a little silly to keep calling it a volume 
so Vj = A Xj. This results in the semi-discrete form of (3.1.1). 

= 0 (3-1.6) 

This equation is exact, as written. Using it requires computing and fj-±- 

Assuming this is possible, the result of a time integration is a collection of values qj at 
some future time. Usually what is needed is q(x) instead, so some assumption needs to be 
made about the distribution of q within each volume. 

Consider integration of (3.1.6) over a time sufficiently small that qj is essentially 
constant. Then it is important to be sure that reconstruction of q(x) from qj is done in 
such a way that entropy is not destroyed. From this point of view the safest assumption 
is the piecewise constant distribution of Godunov. 

q{x) = qj Xj_i <x< x j+ x (3.1.7) 

As shown by (2.2.6), this distribution provides the maximum entropy in each cell, consistent 
with the known values of qj. 

Starting with (2.1.7), a semi-discrete version of the second law can be constructed as 
well. This is 

(P.)i = (S,)„ + > 0 (3.1.8) 

The cell average entropy density, is defined 

S iS l- f S(,)dv (3.1.9) 

v i JVj 
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(3.1.10) 


Because of the assumed piecewise constant distribution of q, 

Sj = S(q t ) = S( ti ) 

and the chain rule is applicable. In particular 

(&it)j = (3.1.11) 

Substituting this into (3.1.8), using (3.1.6) to eliminate (g,t)> - , and multiplying through by 

A Xj gives 

+ (Fj+l < 3 - 112 ) 

This semi-discrete form of the second law has no time derivatives and can be computed 
from data at a single time level. The question of computing such quantities as fj+§ and 
■fj+l must now be addressed. Reviewing their origin (3.1.4), it follows that 

f Hi =f(q i+i ) (3.1.13) 

F j+i =F(q j+i ) (3.1.14) 

where q J+ 1 is the state which exists on the boundary between cells j and j + 1. In virtually 

all of the work which follows, is approximated from qj and qj+i. Definitions (3.1.13) 
and (3.1.14) are then used to compute and Fj + ^. 

The definition of fj+^ and Fj + ^ from a common value of ^ , as in (3.1.13) and 
(3.1.14), differs from the usual practice. Usually / J+ | is formed directly from f(qj) and 
f(qj+ 1 ) without ever computing 1 . The quantity Fj + x is approximated in some similar 
fashion. The author feels that the independent approximation of Fj + i and fj+± does 
not do justice to the second law. Such an approach generates values of which can 
be grossly inaccurate in the neighborhood of shocks, where F is discontinuous. This can 
result in local anomalies such as oscillations. However, because of the global cancellation 
of entropy fluxes, stability is not affected by such errors. Support for this point of view is 
provided in Chapter 9, with a discussion of Tadmor’s scheme. 

It may seem that the piecewise constant assumption (3.1.7) precludes a unique value 
of q J+ i . This can be resolved by relaxing it slightly, in the immediate vicinity of x J+ i , to 
make q(x) continuous. The modified assumption would be 

?(*) = Ij x j-$ +c < « < - « (3.1.15) 

where e « A Xj is small enough that the integrals are not affected and the chain rule still 
holds. 

Because q(x) is continuous, any analysis valid for the integral equations is also valid 
for the semi-discrete equations. Due to the piecewise constant assumption, the domain 
integrals can easily be carried out cell-by-cell as sums. The conjecture is repeated here in 
the semi-discrete form. 
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Conjecture. 


Given 


fj+i - f\ 

(3.1.16) 

J 

(9j)>t + ij+\ - f\ = 0 
i=i 

(3.1.17) 

U)H 

II 

(3.1.18) 

J 

Y ~ ^ ® 

i=i 

(3.1.19) 

_ E /=! 

E 

(3.1.20) 

There exists a value q* which can be determined from q such that 


0<X>- «i) T s4,te - 5/) < E - %>) 

(3.1.21) 


i=i i = i 


The terms having to do with initial conditions have been omitted here to avoid clut- 
tering the notation. Otherwise there is no problem with including them. The constraint 
(3.1.19) can be satisfied with the cellwise condition 

(P.)i > 0 (3.1.22) 


where (P s )j was defined in (3.1.12). 

The major change, in going from the integral form to the semi-discrete form of the 
equations, is the piecewise constant assumption (3.1.7). This means that g(x) is con- 
strained to be piecewise constant for all time. This does not affect the conservation laws 
or the second law of thermodynamics. Stability does not seem to be affected either, since 
the stability conjecture of the previous chapter continues to hold. There is, however, one 
major area where this constraint plays a role. That area is accuracy. 
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3.2 Accuracy in the Fine Mesh Limit 

So far, the claim has been made that schemes which are conservative and satisfy the 
second law will be nonlinearly stable. In spite of this desirable property, such schemes are 
not useful unless they converge to the correct solution. One scheme which does not, is 

(3.2.1) 

This scheme conserves everything and qj°° = qj°. Unfortunately, qj°° is generally incorrect. 

At a mi ni mum , schemes must be consistent in the sense of Lax if they are to converge 
to the correct solution. That is, if a numerical flux function, h, is defined such that 


f}+\ = Mft-i > » Qj+i > 9j+2 ) 

(3.2.2) 

then the consistency condition is 


M 9»9>9*9) = /(?) 

(3.2.3) 


i.e., if the solution is a constant, the computed flux should have no error. This amounts 
to accuracy in the lowest order Taylor series sense. 

Consider schemes of the form 


/ )+} =/( ?j+1 ) (3.2.4) 

where 

9j+ \ = «i+ |(9i» 9i+» ) (3.2.5) 

Consistency for such schemes immediately follows if 

?i+|(9,9) = 9 (3-2.6) 

This condition will be used in most of the work that follows. 

The question of accuracy remains. Traditionally, this question has been approached 
from the point of view of a Taylor series. To simplify the analysis, consider the scalar case. 
For example, (3.2.6) becomes 

«,+$(«* «) = « (3.2.7) 

This work uses interpolations of the form 

U j+ 1 = Uj + (uj+1 - «j) (3.2.8) 
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where <f>j + t is an adjustable parameter, typically 0 < < 2. Interpolations such as 

(3.2.8) are used to form derivative approximations such as 


Ax(u,*) y « U H | - 1 (3.2.9) 
The expression for u j-i is a simple index shift of (3.2.8). Substituting into (3.2.9) gives 


Az(u, t )j « 


uj + 2^;+i(«j+i -«;) 
U i-1 + 


(3.2.10) 


Assuming a constant Ax, and expanding the terms Uj+i and uj-j about Xj gives the 
formula 


Ax 


= («>*); |l + j (^i+i “ ^i-$) 
(«,**)i [l - 

(«,**« )i 1 + ^(<A,+ * ~ «^-a)J 

+0(Ax 3 ) 


Ax 

2~ 

Ax 2 


(3.2.11) 


At first glance, consistency appears to require that <{>j+i = 0j_i. By induction, this 
implies a constant value of <j>. By the same line of reasoning, second-order accuracy seems 
to require that the second term vanish, which happens if — 2. Thus <f> = 1, 

the symmetric case, appears to be the only choice. 


Fortunately, appearances can be deceiving. Less restrictive conditions for second-order 
accuracy can be written 


*,•+}- ^-j=0(A* J ) (3.2.12) 

1 ~ + ) = 0(& x ) (3.2.13) 

This point of view allows a much broader class of second-order accurate schemes to be 
considered. For example, consider the choice 


It may be verified that 


j. - u i ~ u J-i 

,+ > Tj ~ ttj+l “ Uj 

(3.2.14) 

* i+i = r, =1 -xfT i ^ + 0(A*’) 

\ U ^)j 

(3.2.15) 

*j-t = r j—i = 1 - + 0(A*») 

1 (U,z)j 

(3.2.16) 
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although rj Tj- 1 . By inspection, this scheme meets conditions (3.2.12) and (3.2.13) and 

is therefore second-order accurate. This may be independently confirmed by substitution 
of (3.2.14) into (3.2.10). Dividing through by Ax yields the well-known approximation 




Ax Ax 

= («,*)j + 0(A* 2 ) 


3 o , 1 1 

-Uj - 2«j_i + 


(3.2.17) 


Conditions (3.2.12) and (3.2.13) allow some freedom in the selection of <t>j+i and <f>^_ i . 
This freedom will be exploited in succeeding chapters to produce second-order accurate 
schemes which satisfy the second law. 

3.3 Optimal Accuracy — Nonlinear Programming 

Any time the term optimal is used, it must be carefully defined. Section 3.1 defines a 
cell entropy inequality (3.1.22) which is sufficient to guarantee nonlinear stability, at least 
for scalars. Equation (3.2.6) identifies a further constraint (consistency) which is required 
for convergence to the correct solution. Finally, equation (3.2.8) provides the freedom to 
achieve second-order accuracy in spite of these constraints. An optimally accurate scheme 
is defined here as the one which 

(1) has the smallest possible entropy production rate and 

(2) satisfies (3.1.22) and (3.2.6). 

When q is a scalar, the selection of n independent values (q J+ ^) is sufficient to de- 
termine a scheme. For optimal accuracy, these n unknowns must minimize a nonlinear 
functional while satisfying 2n — 1 coupled, nonlinear constraints. This is known as a non- 
linear programming problem. Such problems are thought to be NP complete which is to 
say, they are completely intractable for any reasonable value of n. 

Due to recent developments, the linear programming problem (linear functionals and 
constraints) is no longer NP complete. Solutions for n up to 1000 are now feasible. This 
is still not sufficient for our purpose in more than one dimension, even if the nonlinearity 
could be accommodated. 

To make the problem more tractable, (while waiting for further advances in solution 
of the nonlinear programming problem) strict optimality can be given up in exchange for 
uncoupling the constraints. The following section demonstrates how this might be done. 

3.4 First-Order Accurate Schemes 

The first challenge in this nonlinear optimization problem is to show that the con- 
straints can be met. This section will prove that they always can be, at least in the scalar 
case. 
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One constraint is the semi-discrete cell entropy inequality. This constraint is nonlinear 
and is applied at each control volume. Consider here the semi-discrete form of the second 
law in one dimension. 

(.Si),, + > 0 (3.4.1) 

This was originally derived in Section 3.1. There it was also shown (3.1.12) that an 
equivalent statement is 


+ (3.4.2) 


Equation (3.4.2) is algebraically identical to 

[-($.,),(/,•+ j - Si) + (Fi+i - */)] + [-(S..M/, -fi-i) + (Fi - F M )] > o (3.4.3) 

The left-hand side is an expression for the total rate of entropy production in the cell, 
(P,)j. The terms in brackets define the right and left half cell entropy production rates, 
(P+ )j and (P”)j. The dependent variables here are fy+j (which defines / J+ 1 and 
and 9j-i (which defines fj-x and F J+ ^). As written, (P?)j depends only on 1 , and 
(P“)j depends only on qj_ i. This suggests replacing (3.4.3) with the more restrictive 
constraints 

(P+)j >0 and (p-)i > 0 (3.4.4) 

Adoption of (3.4.4) allows Qj+x to be selected independently of This decoupling of 

the constraints makes solution feasible. 

Constraints (3.4.4) imply that Qj+x must be chosen to simultaneously satisfy the 
following two inequalities 


(P, + )i> 0 and (P,-)j+i > 0 (3.4.5) 

since (3.4.4) must be satisfied both at point j and point j + 1. Is this always possible? For 
scalars, the answer is yes. For systems, the answer is not known in general. The argument 
follows. 


Begin with an initial guess for a simple average of qj+i and qj. This allows 

(P+)j and (P~)j+i to be evaluated. If either of these is negative (the usual case), then 
qj+x will need to be adjusted. Say that (P+)j is negative. Then an appropriate change to 

qj+x is along the gradient, . The differentiation is easily carried out 


djPTh 


— (SMf w)j+i 


(3.4.6) 
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(3.4.7) 


Using the compatibility condition (F„ = S, q f, q ) this can be rewritten 

^2]* = {(S.,)i+) - (5„)>}(/.,) j+ i 

Adjusting <}j+± in this direction will increase (P+ )j as planned, but the same adjustment 
may decrease (P~)j +i . This can be prevented if is restricted to values for which 


( s(p. + )A \ ^ 

V / V ) 


0 


(3.4.8) 


It is instructive to carry out the differentiation and the dot product using (3.4.7). 

(K in) ( d(p.-) i+ , \ 

V 8 «i+i } \ dq i+i ) (3.4.9) 

Using Taylor’s theorem with remainder 

(*«)/+* ~ {S n )i = {q j+h - q,)S- q (3.4.10) 

(for some q m ) which allows (3.4.9) to be written as 

( ' (^r) = Ui+ i ~ - W (^U) 

If q is a scalar quantity, then S y * q and 5,** are negative scalars (from convexity). 
Furthermore, (f,q)j+x is a scalar. It follows that condition (3.4.8) is satisfied for scalar q, 
if and only if 

(qj+i - q i+ \ )(«;+! - qj) > o (3.4.12) 

In words, this says that must lie between qj and qj+i- Within this range, any 
adjustment to q J+ i which increases (Pf)j will not decrease (P~)j+i. 

As long as does not change sign, can be adjusted all the way to qj at 

which value (P+ )j = 0 by inspection. If (P“)j+i is negative, qj+x can be adjusted to 
qj+ 1 . In either case (3.4.1) will be satisfied. 

Under what conditions can (and 0 ) change sign? Substituting 3.4.12 

into 3.4.9 gives 


= (hh - 


(3.4.13) 


27 



Because of convexity, S,* q is a negative constant. Because of (3.4.12), Qj+i — qj has the 
same sign as qj+i—qj. Only (/, 9 ) J+ i can change sign on the interval. If this happens, there 
will be a value of qj+i on the interval, at which (/,,) J+ 1 = 0. If this represents a local 
maximum (sonic point), both (P+)j and (P^)j+i will be positive there. If it represents a 
local minimum (shock), will be adjusted away from it, to the appropriate endpoint 
value. 

As a kind of side effect, the consistency constraints have been satisfied. Notice that 
(3.4.12) is sufficient to satisfy (3.2.6), which itself is sufficient for consistency in the sense 
of Lax. 

A typical first-order upwind scheme might use the formula 

q j+ ! = qj + ^(1 - sign((/„)j + (/„)j+i))(g;+i - qj) (3.4.14) 

everywhere. The only change proposed here is the treatment of sonic points. At these 
points the formula 

— q»o nic (3.4.15) 

is used, where q IO nic is the value at which (P+)j (and (P~)j+i ) reaches a maximum. 

The scheme described by (3.4.14) and (3.4.15) is first-order accurate in the limit of 
mesh refinement. The value <f> referred to in Section 3.2 will always be 0 or 2 (except at 
sonic points), independent of the value of r. To achieve second-order accuracy, the scheme 
must be modified in such a way that <f> — ► 1 as r — *• 1. The next section will show how this 
can be done. 


3.5 Second-Order Accuracy 

In the previous section, an expression for the semi-discrete cell entropy inequality was 
derived. This is restated 

(ft); EE + (P+)j > 0 (3.5.1) 

The quantities (P”)j an d (P, + )j were defined shortly after (3.4.3). In Section 3.4 the 
constraint (3.5.1) was satisfied term by term. This is not necessary. 

Begin as before by choosing, as an initial guess for qj+i , a simple average of qj and 
qj+ 1. When this is done, trial values of (P~)j and (P+ )j can be computed for each control 
volume. For illustration, suppose that (P“)j = 10and(P+)j = —9. The first-order scheme 
in Section 3.4 would adjust qj+i until = qj , in order to assure that (P+ )j > 0. The 
resulting entropy production for the cell, (P«)j, would then be at least 10. In fact, no 
adjustment is necessary. Using the trial values, (P t )j = 1 which meets (3.5.1). 
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On the other hand, suppose that (P+)j = —11 instead of —9. This would have no 
effect on the first-order scheme of the previous section. The second order scheme would 
find it sufficient to adjust Qj+i very slightly, so that (P+ )j = —10 instead of the trial value 
of —11. With this change, (Pg)j = 0, which satisfies (3.5.1). 

In both of these cases, satisfying (3.5.1) as a whole instead of term by term results 
in an order of magnitude reduction of the cell entropy production rate (dissipation). As 
the grid is refined, the trial values of (P~)j and (P+ )j come closer and closer to satisfying 
(3.5.1). The needed adjustments get smaller and smaller. In the limit as the mesh is 
refined, the adjustments become vanishingly small, and —> \(qj + qj+x). In this limit 
the scheme becomes second-order accurate. 

Finally, suppose that (P+)j = —10 and (P~)j = 0. In this case, must be 
adjusted until (P+ ) ; - = 0. In this case (which generally occurs at extrema) the large 
change in cannot be avoided. The scheme reverts to first-order accuracy at such 
points. Generally the number of extrema in the solution is small and does not increase 
with the number of mesh points. 

A convincing demonstration of second-order accuracy requires a Taylor series expan- 
sion (without remainder this time) of (P,)j about the point j. This immediately leads 
(after dropping the constant terms) to 

(P»)j w (P>g “ &iqf >g)(?j+$ — Qj) + {F } q — Syqf jj)(9j_l — Qj) 

+^(9;+! - Qjf( F m - S iq fyqq)(q j+i - qj ) (3.5.2) 

— — — S>qfiqq)(Qj-% — Qj)> 0 

where everything is evaluated at point j unless otherwise noted. The first two terms vanish 
because of the compatability condition, Fyq — Syqfiq. Using this relation allows the second 
two terms to be simplified by expanding F m . The simpler form is 

(P#)j * (Qj+$ — Qj) (Smf*i)(Qj+$ ~ Qj) ~ (?i-| — Qj) (Swf>q)(.Qj-$ ~ Qj) — ® (3.5.3) 

Using (3.2.8), the unknown values of Qj+i can be replaced with unknown values of 4>j+i 
and After multiplying by 4 this gives 

(P«)i ~ Qj) (Syqqf >q)4 > j+ — Qj)~ 

( Qj-i ~ Qj) (&>qqf’q)(2 ~ &j-$) 2 (Qj-i ~ Qj) ^ 0 

Using (3.2.11), the quantity (qj-i — Qj) can be replaced with rj(qj+x — qj). 
terms (valid for scalars) gives 

(P,)j W ( Qj+l ~ Qj) S, qq f,q{<l>j+^ — Tj 2 ( 2 — <j>j_ l ) 2 ) > 0 


(3.5.4) 
Collecting 

(3.5.5) 
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The initial guesses for correspond to <f>j+i = <f>j~ 1 = 1. As the mesh 

is refined, rj — ► 1 . As this limit is approached, the values of <f> will require less and less 
correction from the initial guess. Their final values will approach one. This is the condition 
given in Section 3.2 for second-order accuracy in the limit of mesh refinement. 

In most solutions, there will be some regions where the second derivative cannot be 
neglected and r will not approach 1, no matter how much the mesh is refined. Typically 
these include discontinuities and, in some cases, extrema. At such places, a Taylor series 
shows first-order accuracy. 

The second-order scheme differs from the first-order scheme in that (3.5.1) is satisfied 
as a whole instead of term by term. Shocks and sonic points are not affected. For the 
present, the scalar case can be considered complete. Chapters 5 and 6 give examples and 
results for some common scalar test problems. In Chapter 7, the formal extension to the 
Euler equations will be given. 
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Chapter 4. THE EFFECT OF TIME ADVANCE SCHEMES 
ON ENTROPY PRODUCTION RATES 


In discussing the effects of a finite time step it may help to recall a similar problem 
for the incompressible Navier-Stokes equations. When these were first being solved, there 
was considerable difficulty with stability. Using spatial difference schemes with very high 
order spatial accuracy does not help. One solution to the stability dilemma is to con- 
struct schemes which satisfy a global inequality for kinetic energy. Except for the effect 
of boundary conditions, kinetic energy decreases steadily. Experience with incompressible 
Navier-Stokes solvers is instructive because kinetic energy is almost a formal entropy for 
these equations. Generally, such analysis is semi-discrete; i.e., it is done under the assump- 
tion that the time advance is analytic. This has seemed to work fairly well, even though 
time advance is generally not analytic. 

Now the difference between a semi-discrete scheme (continuous in time) and certain 
fully discrete schemes (discrete in time) can be formally analyzed. In particular, the 
entropy production rates can be bounded on a cell-by-cell basis for both explicit Euler and 
implicit Euler time advance. For scalar equations, a modified form of Crank-Nicolson can 
also be analyzed. 


4.1 The Fully Discrete Case 

In the previous chapter, the integral equations were given on a volume-by-volume 


basis. 


/ q[t + At]dv — f g[<]du+ f <£ f-ndAdr = 0 (4.1.1) 

JVj JVj Jt JdVj 


To simplify notation, time is indexed by a superscript so that 


9? = 9j[t ] 9? +1 = <lj[i + At] (4.1.2) 

Using the definitions (3.1.3) and (3.1.4), (4.1.1) becomes 

/ t+At ( f x — f . t \ 

A*/ ' dr = » (4-1.3) 

which is the fully discrete analog of (3.1.6). Like (3.1.6) it is exact. All that remains is 
to integrate the fluxes over the time interval. Since the entropy flux will also have to be 
integrated, a simple strategy is to assume that fy+i remains constant over the entire time 
interval. For example, explicit Euler time advance is equivalent to 


for all nAt < r < (n + 1)A t 


q H i[r] = q^x 
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(4.1.4) 



and implicit Euler is equivalent to 


qj + 1 [r] = for all nA< < r < (ra + 1)A t 


(4.1.5) 


With the addition of definitions like 


f?+i*W + »> 


the fully discrete form of the equations becomes, for explicit Euler 


(4.1.6) 




(4.1.7) 


Using a similar derivation, the fully discrete form of the second law becomes 




(4.1.8) 


Equations (4.1.7) and (4.1.8) will be used in the next section to show that schemes 
using explicit Euler time advance have a lower (perhaps negative) entropy production rate 
than the same scheme with analytic time advance. This is true, cell by cell, for any positive 
value of At. 

4.2 Explicit Euler Time Advance 

This section explores the irreversibility of explicit Euler time advance. The entropy 
production due to time advance can be separated from the entropy production due to space 
differencing. This is shown in the following theorem. 

Theorem. 


Given a fully discrete scheme using explicit Euler time advance 

„ /;+i - 

(«r - 4?) + At , , Ax . 1 = 0 


Then 


where 


(P.) f* 5 


( P.)f* < (P.h 

_ 5 "+> - 57 r; +i -r;.i 


( 4 . 2 . 2 ) 


( 4 . 2 . 3 ) 
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and 


( P.)i = (S„)j + 


pn _ pn 
r j+l r j~ 1 
J ^3 l 2. 


A Xj 


fn __ fn 
/-•+i /-•_! 




Axj 


(4.2.4) 


Proof: 

Even though (4.2.4) is semi-discrete, q(x) must be specified. Here, q(x) at time level 
n has been selected. The quantities 5” +1 and 5J 1 are related through Taylor’s theorem 
with remainder: 


s? +1 = 5? + (S„)?(?; +I - *?) + - ?7) (4.2.5) 


where n < £ < n + 1. Equation (4.2.5) is exact, however the value of £ is not known. The 
right-hand side of (4.2.5) may then be substituted for 5" +1 in (4.2.3) to give 


(A)?* = (( S ? + ( S „)?(«? +1 - «?) + 5(4” +1 - 1 ? f S 4(4" +1 - «*)) - S ?) 


Ax, 


Furthermore, equation (4.2.1) can be used for the value of (q” +1 — qj). This gives 

{ fn _ fn pn pn ^ 

+ ■ i+ L, H } 


(4.2.6) 


+ 2S7 ( «” + ’ - «") Ts ’«(«i +1 - «?) 


(4.2.7) 


The quantity in braces can be replaced with the semi-discrete entropy production rate 
through (4.2.4) to yield 


(P,)f E = (P,)j + ~ ~ *’) (4-2.8) 

Since £,35 is a negative definite matrix (by convexity), the second term of (4.2.8) is a 
proper quadratic form and cannot be positive. Thus, the fully discrete entropy production 
rate (P t )f E must be equal to or less than the semi-discrete entropy production (P a )j- 

What does this say about explicit Euler as a time advance scheme? For problems 
which naturally create entropy (such as the heat equation), it may be a reasonable scheme. 
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Since (qj +1 — g”) is roughly proportional to At, the quadratic term is first-order in time. 
This agrees with existing analysis that explicit Euler is first-order accurate. 

Stability depends on the relative size of the (presumably positive) semi-discrete en- 
tropy production rate, (P a )j, and the (generally negative) quadratic term. For small 
enough time steps, (P s )j dominates, and the fully discrete entropy production rate is pos- 
itive (( P,)f E > 0). For very large time steps the opposite occurs. Thus the entropy 
approach demonstrates that explicit Euler is conditionally stable for dissipative systems. 
As a steady state is approached', (q” +1 — qf) vanishes and the entropy production rate is 
entirely determined by the space differencing, as expected. 

Convection problems are the main concern in this work. Provided that the derivatives 
exist, such problems can be expressed as 

9,t + /,E = 0 (4.2.9) 

If q(x,t) is smooth enough to allow use of the chain rule, it can easily be shown that 
P a = 0, analytically. This was shown in Section 2.2. 

The more accurate a space differencing scheme becomes, the closer it comes to this 
analytic limit. An ideal space differencing scheme, one with no error, has no semi-discrete 
entropy production, except at shocks. For such a scheme, explicit Euler is not stable for any 
positive time step At. The usual fix is to add artificial spatial dissipation, which modifies 
the equations. The resulting solutions can be surprisingly good for transient problems, if 
the spatial dissipation (an error) approximately cancels the temporal dissipation (an error 
of opposite sign). The first-order upwind scheme, for example, produces no error for the 
wave equation if At = cAx (where c is the local wave speed). This strategy of canceling 
errors is difficult in practice (for all but the simplest problems) because c is not constant in 
space or time and because there may be more than one characteristic speed. Furthermore, 
it doesn’t work at steady state or at very small time steps because the time advance errors 
vanish, while the spatial errors which they are supposed to cancel, remain. 

4.3 Implicit Euler Time Advance 

This section explores the irreversibility of implicit Euler time advance. The entropy 
production due to time advance can be separated from the entropy production due to space 
differencing, just as for explicit Euler. This is shown in the following theorem: 

Theorem. 

Given some fully discrete scheme using implicit Euler time advance 

f n + 1 _ yn+1 

(«" +1 - 1?) + AiJ±i ^ t± = 0 (4.3.1) 
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Then 

{P.)] E > (P.)j 

(4.3.2) 

where 

Cn+l _ cn P n + 1 - F n+ ? 

1 t)j ~ A t 1 Ax,- 

(4.3.3) 

and 

pn+t pn+1 

(P.)i = (5..)? +1 + 



/n+1 /n+1 pn+1 rm+1 

(4.3.4) 


Proof: 


Even though (4.3.4) is semi-discrete, q(x) must be specified. Here, ^(x) at time level 
n + 1 has been selected. The quantities S” +1 and Sf are related through Taylor’s theorem 
with remainder: 

5 ; = s? +i - - «”) + 5<«" +1 - «") Ts .;,(«? +1 - «?> («•*) 

where n < rj < n + 1. Equation (4.3.5) is exact, however the value of rj is not known. The 
right-hand side of (4.3.5) may now be substituted for 5” +1 in (4.3.3) to give 


(p.Yj E = (s? +l - (s ? +1 - (s„) 7 +, (<? +1 - q?) 

jpn+l _ j?n+l 

+ (4.3.6) 

Furthermore, equation (4.3.1) can be used for the value of (q? +1 — qf )- This gives 


(P.)', E = < 


/n+1 /n+1 cm+1 — pn+l 1 

. ^ J — "i — — i+1 j — — I 

— {S 2 J 2 _|_ ~2 J _2_ Vu 


Axj 


Ax. 


+ ^(47 +I -4" )T S.:, ( 4? +I -S" ) 


(4.3.7) 


The quantity in braces can be replaced with the semi-discrete entropy production rate 
through (4.3.4) to yield 


(p.)i B = (A),- - - 9?) 


(4.3.8) 
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Since S, qq is a negative definite matrix (by convexity), the second term of 4.3.9 is a proper 
quadratic form. Thus, the fully discrete entropy production rate ( P,)j E must be equal to 
or greater than the semi-discrete entropy production ( P $ )j. 

What does this say about implicit Euler as a time advance scheme? Since (?” +1 — qf) 
is roughly proportional to At, the quadratic term is first-order in time. This agrees with 
existing analysis that implicit Euler is first-order accurate. In this regard implicit Euler is 
very similar to explicit Euler. 

In the area of stability, though, the two schemes are very different. Stability depends 
only on constructing a space differencing scheme for which (P t )j > 0. It then follows that 
(P»)j E > 0 also, independent of the time step. Thus the entropy approach demonstrates 
that implicit Euler is unconditionally stable in a nonlinear sense for convective or dissipa- 
tive systems. The catch is that this wonderful property only follows if the scheme is carried 
out exactly. For most problems, this implies solving coupled systems of nonlinear equa- 
tions, an expensive process. Solution of these systems is possible if the time step is small 
enough to allow linearization in q, however. Thus, while implicit Euler offers unconditional 
stability, this advantage is fully realized only in linear problems. For nonlinear problems, 
some time step restriction may result from linearization errors. This restriction relaxes as 
a steady solution is approached. Furthermore, at steady state, (qj +1 — qj) vanishes and 
the entropy production rate is entirely determined by the space differencing, as expected. 

4.4 A Second-Order Accurate Time Advance Scheme 

It is interesting to note the similarity in the analysis of explicit Euler and implicit 
Euler time advance. Implicit Euler is essentially the same as explicit Euler with a negative 
time step. The errors incurred (comparing (4.3.8) with (4.2.8)) are of the same form, but 
of opposite sign. This suggests that it might be possible, by combining the two methods, 
to construct a time advance having no temporal dissipation. In particular, if a half step 
of implicit Euler is followed by a step of explicit Euler, the temporal dissipation might be 
partially or completely eliminated. 

Theorem. 

Given some fully discrete scheme; a time advance of the form 

f n +\ _ f n +\ 

(«? +1 - 1?) + A t = 0 (4.4.1) 

which is composed of a half step of implicit Euler ( t <t <t + followed by a half step 
of explicit Euler (t + ^ < r < t + At ), 
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Then 


(P.)? ON 


= (^ + 8 A ?« 


r»+l 




(4.4.2) 


where the precise evaluation states of S,% q and S,* g are not known. Equation (4.4.2) uses 
the definitions 


(P.)f CN = 


S"+i _ S v. 

At 


F n+ > - F n+ ’ 

j+i jzi 

Axj 


(4.4.3) 


and 


F n +i _ ^n+i 

(p.)i s (S„)" + ’ + 


3~1 


Axj 




f n +5 _ f n +l 

; _M 1M 


A xj 


+ 


F-*} - F" + , } 

3_+j £zi 

Axj 


(4.4.4) 


Proof: 


Here, q(x) at time level n + | has been selected for the semi-discrete terms. The 
quantities and 5” are related, through Taylor’s theorem with remainder to : 


s? = - ( 5 „>r’ <«; + ’ - *?) + ^r 4 - 1 - «?> 


(4.4.5) 


s; +1 = «r* + (s„)” +1 (,; +i - ,;+*) + i(,; +i - «” +i ) T s, s,(«? +1 - «; +i > 


(4.4.6) 


where n<7;<n + |<£<n + l. Substituting (4.4.5) and (4.4.6) into (4.4.3) gives 


(P.)j 


MCN _ 


+ 5(?? +i - ?r , ) 7s 4(«r i - «r +i )> 

-(s; + i -(s„); + !(,; + i 

+ -«;») 


+ 


p n+ j _ F n+ i 

i+i 


i 2 


Axi 


(4.4.7) 
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Collecting terms gives 


(p.)Y cN = ^((s„)r*(»r' - <n + - «; +, )\ £ ,(«" + ‘ - «; +i ) 

- - i?) T s,iM +i - «”>) 


+ 


- F n+ * 

j+i * - 1 


J : 


A Xj 


(4.4.8) 


The key observation at this point is that 


<9? +1 - 9,” +i ) = (jf’ - 9?) = -f /j>! Az 7' L = " «"> 

Furthermore, equation (4.4.1) can be used for the value of (q” +1 — 9™)* This gives 


n+ 


(4.4.9) 


(A)f CN = < -(5„); 


n+$ J j+\ , J+j M 


Axw 


+ 


Ax j 


+ 8^(9" + ' " 9"f (S4 - «.;,)(9r ' - «“) < 4 - 41 °) 

The quantity in braces can be replaced with the semi-discrete entropy production rate 
through (4.4.4) to yield 

(P.)f ' ™ = (P.)i + 8^(9" +1 - 9 7f( s 4, ~ S,;,)(9” +1 - 9?) (4-4.11) 

which is the claimed result. 


Although S,t q and S,% q are both negative definite, the sign of their difference is inde- 
terminate. Although they are evaluated at states which may be different, both are, in a 
Taylor series sense, evaluated close to • Their difference is of order At. This accounts 
for the second-order accuracy in time that this scheme enjoys under conventional analysis 
techniques. 

The behavior of this time advance scheme for scalars is quite remarkable. As previ- 
ously noted, S,t q = S^ q = —2 in this case. Equation (4.4.11) indicates that ( P t )j^ CN = 
(P»)j ; the scheme provides no temporal dissipation at all. This is true cell-by-cell. Thus, 
the discrete time advance is as good as an analytic time advance. This follows regardless 
of the time step. 

Since half the scheme is implicit Euler, it is reasonable to expect that the linearization 
problems of implicit Euler will carry over. Indeed they do. The scheme is, however, no 
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more difficult to implement than implicit Euler. Given that one can solve for qj using 
implicit Euler, the equivalent higher order scheme is 


?" +1 - 


= + 2(?7 +i - 9 ?) 


(4.4.12) 


This follows directly from (4.4.9). Many codes in use today use the so-called “delta” 
form. Equation (4.4.12) shows that, simply by doubling the correction (a sort of overre- 
laxation), it is possible to achieve second-order accuracy in time. This fact is apparently 
known, but seldom used. 

There is a connection between this scheme and the more common Crank-Nicolson, 
time advance. Crank-Nicolson consists of a half step of explicit Euler, followed by a half 
step of implicit Euler. Though the scheme shown here reverses the two, a sequence of 
Crank-Nicolson steps preceded by a half step of implicit Euler and succeeded by a half 
step of explicit Euler will yield exactly the same result as the scheme suggested here. 

Traditionally one includes, for generality, the so-called 9 schemes in an analysis of this 
sort. These are variants where one advances in time by 9 At, using implicit Euler and then 
advances (1 — 9) At, using explicit Euler. Values of 9 = 0, and 1 correspond to explicit 
Euler, modified Crank-Nicolson, and implicit Euler. These are the only values ever used 
in practice. Suffice to say that the analysis shown in this section is easily extendable to 
other values of 9. 

4.5 Practical Considerations for Implicit Schemes 


In linear problems, it is possible to implement implicit schemes with style and grace, 
retaining rigor and exactitude. On the other hand, most problems of interest are nonlinear. 
What can be done in these cases? Usually all that can be done is to locally linearize and 
take small time steps. This section will address the question of how to maintain the entropy 
inequality while doing this. 


Implementing implicit Euler time advance (eq. 4.3.1) requires evaluation of terms 
such as fj+i • In general, this is very difficult. What is easy to evaluate is If the 

two fluxes are not too different, it is reasonable to use a Taylor series in place of the 
unknown quantity. Supposing, for illustration, that = f{qji9j+ 1)> then a natural 

approximation is 


fn+l 

f j+\ 



+ 


dqj 




~ 9?+i) 


(4.5.1) 


This approximation makes implicit methods possible. It reduces a coupled set of nonlinear 
algebraic equations to a coupled set of linear algebraic equations. Such systems are easy to 
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solve, expecially in one dimension. Typically a tridiagonal or pentadiagonal matrix needs 
to be solved. 

How appropriate is this approximation? If (q” +1 — qf) and {qj+i — 9J+1) are “small 
enough,” then it is an excellent approximation. Their size is, of course, directly related 
to the time step At. By controlling A t we can make them “small enough.” How small 
is that? Small enough that the discrete entropy production rate (e.g., eq. 4.3.8) remains 
positive. 

Another problem appears. Even for a very small (finite) time step, there are lineariza- 
tion errors in the approximation (4.5.1). This introduces errors, even in an a posteriori 
calculation of the entropy production rate. These uncertainties jeopardize the conclusions 
of Section 4.3 (and Section 2.5). 

The most straightforward solution is to iterate within each time step to exactly solve 
the nonlinear algebraic equations. This is expensive however and may not be robust due 
to certain roundoff difficulties. 

A cheaper plan is an implicit-explicit scheme of the following type. First perform an 
approximate implicit step to get an approximation to q n+1 . 

A (q~q n ) = S x f (4.5.2) 

where A is the matrix implicit in (4.5.1) and S z f is a vector of flux differences. The vector 
q is approximately (exactly for linear problems) equal to q n+1 . Equation (4.5.2) allows 
solution for q at all grid points by standard linear algebra techniques. 

Second, compute the fluxes for this approximate solution. 

(4.5.3) 

The quantities qj+i are computed from qj and q J+ i using whatever semi-discrete scheme 
is chosen. 

Finally, use these fluxes explicitly to compute q n+1 . 

This allows an a posteriori check on the entropy inequality everywhere. 

%? +1 ) - *(«?) + (^W - })) ^ 0 < 4 - 5 - 5 ) 

If it fails anywhere, the time step is too large. With this test, one can be sure that ap- 
proximation (4.5.1) is “good enough,” at least in the sense that stability will be preserved. 
A similar strategy can be performed for the second-order accurate scheme of Section 4.4 
through the implementation shown in equation (4.4.12). 
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Chapter 5. LINEAR SCALAR EQUATIONS 


5.1 A Suitable Entropy for the Wave Equation 

This chapter will consider solution of the ordinary scalar wave equation in one di- 
mension. Despite its simplicity, this equation provides a good illustration of the principles 
involved. It is a reasonable model for a contact discontinuity, because in both cases the 
characteristics are parallel. 

For this and other scalar equations the variable u will be used instead of q, mostly 
because these equations are usually written that way. The scalar wave equation in one 
dimension is 

u, t + cUyx = 0 (5.1.1) 

where the constant c is the wave speed. The exact solution is a translation of the initial 
conditions with speed c. The boundary conditions are considered to be periodic. The 
exact solution has no dissipation; it conserves entropy. Numerical schemes satisfying the 
second law will usually produce some (hopefully small) amount of entropy. This results in 
smearing of the initial conditions. 

The results of Chapter 4 leave only the task of finding a spatial differencing scheme 
which satisfies the semi-discrete entropy inequality. Using such a scheme with Modified 
Crank-Nicolson or Implicit Euler time advance produces a scheme satisfying the fully 
discrete entropy inequality. 

Consider first the semi-discrete approximation. 

u, t + = 0 (5.1.2) 

For this equation / = cu. Since c is constant (5.1.2) reduces to 

+ = ° ( 5 . 1 . 3 ) 

The idea here is to pick 1 in such a way as to satisfy a cell entropy inequality. 
Before that can be done an entropy has to be defined. In Section 2.5 it was shown that 
S = —u 2 is a suitable entropy for scalar equations such as this. The convexity condition 
is easily verified. 

5, uu = -2 < 0 (5.1.4) 

A corresponding entropy flux can be obtained by direct integration of the compatibility 
condition 

F iU = S, u f, v (5.1.5) 
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The right-hand side evaluates to — 2cu which is easily integrated to find the entropy flux 
F. A suitable entropy pair then, is 

S = -u 2 and F = -cu 2 (5.1.6) 

This choice allows a semi-discrete entropy inequality to be formulated. Substituting S, u 
for S , q , cu for / and —cu 2 for F, the semi-discrete entropy inequality (3.1.12) is, after 
multiplication by A®, 

A x (P,)j = 2cuj(u j+ i - u ; _i) - c((u i+ i) 2 - (uj-i) 2 ) > 0 (5.1.7) 

Remarkably, this expression factors neatly into a difference of two squares. 

A x (P,)j = -c [(« j+ i - ujf - ( Uj - uj_i ) 2 > 0 (5.1.8) 

This expression for the semi-discrete entropy production rate is exact. All that remains is 
to choose and Uj_ a in such a way that (5.1.8) is satisfied. 

5.2 A First-Order Scheme 

Construction of first-order accurate schemes which satisfy a cell entropy inequality 
was discussed in Section 3.4. The general idea is to eliminate negative terms from the 
entropy production rate. This is easily accomplished in this case. For positive values of c, 
the choice 

Uj+i = (5.2.1) 

will leave only the positive term, c(uj — u j-i ) 2 > in (5.1.8). If c happens to be negative, the 
definition u,_ i = uj should be used instead, leaving the positive term — c(u,- , i — uj) 2 . 
Since c is a constant, never changing sign, there can be no conflicts requiring Uj + i to equal 
both Uj and Uj+i. 

The total entropy production rate for this scheme can easily be calculated. The 
entropy production for the , 7 th cell is 

A x(P t )j = c{uj - Uj_ i) 2 (5.2.2) 

for c > 0. By definition (5.2.1) Uj-i = Uj - 1 in this case, so 

A*(P,)j = c(uj - Uj-i) 2 (5.2.3) 

Thus the total entropy production is the sum of the squares of the differences between 
grid points, a result which holds even on irregular meshes. Since total entropy production 
should be zero, a reasonable measure of error can be obtained by summing A x(P t )j over 
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the mesh. This measure of the total error will always be a positive quantity for nontrivial 
initial conditions. Assume for now that an “optimum” mesh has somehow been chosen, 
one in which each of the terms (5.2.3) is of equal magnitude. This minimizes the total 
entropy produced for a fixed number of cells. The effect of doubling the number of cells 
will be to double the number of such terms, while reducing their individual magnitudes by 
a factor of four. The net effect is to cut the entropy production in half. For this problem, 
the entropy production rate is inversely proportional to the number of grid points, other 
things being equal. This is characteristic of a first-order accurate scheme. 

This is a first-order upwind scheme. It is well known in the literature and is highly 
dissipative unless an entropy destroying time advance scheme is used. For the special case 
of explicit Euler time advance with time step At = ^, the errors in the time advance 
scheme precisely cancel the errors in the space differencing scheme and no entropy is 
produced. At larger values of the time step, the time differencing errors dominate and a 
net destruction of entropy results. The scheme is then unstable. 

5.3 A Second-Order Scheme 

Construction of second-order accurate schemes which satisfy a cell entropy inequal- 
ity was discussed in Section 3.5. In this section the wave equation is used to illustrate 
the technique. The semi-discrete formula for entropy production (5.1.8) is restated for 
convenience. 

Ax(P,)j = c(uj - Uj _ i) 2 - c(u i+ i - Ujf > 0 (5.3.1) 

Assume c > 0. A useful strategy is to reduce the magnitude of the negative (second) term. 
If nothing is known about the magnitude of the first term, then the second term has to 
be reduced to zero to assure satisfaction of (5.3.1). That strategy leads to the first-order 
upwind scheme explored in the previous section. To improve on it, consider schemes in 
which 

I u i - «j-i I > g I u i ~ u i~i I (5.3.2) 

This convention establishes a (usually) nonzero lower bound on the magnitude of the first 
term in (5.3.1). The magnitude of the second term need only be reduced to | \uj — Uj- 1 |. 
Given the constraint (5.3.2), equation (5.3.1) reduces to 

K+i — u j\ < \ Wj ~ Uj-i I (5.3.3) 

The second-order accurate linear approximation 

u j+ i = i(u i+1 + uj) (5.3.4) 

may satisfy (5.3.3). If it doesn’t, a different second-order approximation can be used 
instead 

u j+ i = -uj - -UJ-J (5.3.5) 
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The two formulae are equivalent if (ttj+i — Uj) = (uj — uj- 1 ). This can be summarized as 
follows 


u i+\ 


u j + 2^' +1 “ 


Uj + ~ 


ttj ~ u i-i 
»i+ 1 “ u i 


(«i+i - u i ) 


if |«j+i - tij| < | uj - Uj-i\ 

otherwise 


(5.3.6) 


As in the first-order scheme, no conflicts arise. If — Uj\ is reduced, |«j+i — Wj+i | 

will increase, allowing (5.3.2) to be satisfied at grid point j + 1. 


The commonly used MINMOD limiter differs from (5.3.6) in its use of the first-order 
formula Uj + i = uj at extrema. At extrema, (5.3.6) behaves slightly better than the 
MINMOD limiter. The first reference to (5.3.6) is in an earlier publication by the author 
(ref. 25). 


As in the first-order scheme, the total rate of entropy production can be computed. 
This computation requires understanding the switch given in (5.3.6). If Uj+i = |(uj + 
Uj^-i), the first branch, the scheme is exactly central differencing. The second branch 
evaluates to Uj + i = | Uj — \uj~ \ except at extrema. This corresponds to a second order 
upwind scheme of Warming and Beam. 



Consider the special case where u(x) is the smooth pulse shown in figure 5.1. In 
regions I and III, a second-order upwind scheme is used. In these instances, the quantity 
(uj _ |_ i — Uj), which appears in the entropy production rate, is given by 


(« i+ l - Uj) = 



(5.3.7) 


and the quantity ( Uj 
by 


which also contributes to the entropy production is given 
( Uj — Uj_ i ) = ^(2uj — 3uj_i + Uj- 2 ) (5.3.8) 


44 



When the entropy production is computed according to (5.3.1), using formulas (5.3.7) 
and (5.3.8), the result factors neatly 

Ax(P # )j = ^(3ttj - 4uj_i + Uj- 2 )(uj- 2 - 2uj-i + Uj) (5.3.9) 

The first factor is a first derivative formula (except for the A*) and the second factor is a 
second derivative formula (without the A* 2 ). In fact, a Taylor series shows (dividing out 
the Ax) 

(P,)j = ^u, x u, xx Ax 2 + 0( Ax 3 ) (5.3.10) 


In regions II and IV, a second-order central differencing scheme is used. In these 
regions, the quantity (u J+ 1 — tty), which appears in the entropy production rate, is given 
by 

( u j+\ ~ u j ) = 2 ( u i+t — u i ) (5.3.11) 

and the quantity ( Uj — which also appears is given by 


( Uj -Uj_i) = -(uj-Uj-t) 


(5.3.12) 


When the entropy production is computed according to (5.3.1), using formulas (5.3.11) 
and (5.3.12), the result again factors neatly 


Ax(P,)j = --(«;+ 1 - «i-i)(ttj-i - 2uj -f u j+1 ) 


(5.3.13) 


The first factor is a first derivative formula (except for the A*) and the second factor is a 
second derivative formula (without the Ax 2 ). In fact, a Taylor series expansion shows 

(P.)j = ~^u, x u, xx Ax 2 + 0(Ax 3 ) (5.3.14) 

In equations (5.3.10) and (5.3.14) u, x and u, xx are evaluated at point j, for which the 
entropy production is being calculated. Notice that, in both cases, the entropy production 
rate declines as the square of the mesh size. This implies second-order spatial accuracy. 

Notice that (5.3.10) and (5.3.14) are identical except for sign (and higher order terms). 
This is because the product u, x u, xx is positive in regions I and III while it is negative 
in regions II and IV. Thus the entropy production rate (P t )j is positive in either case. 
Incidentally, equations (5.3.9) and (5.3.13) both evaluate exactly to zero for linear data. 
This is another test for second-order accuracy. 


The remaining cases to consider are the transitions. There are three different kinds. 
Between regions I and II or between regions III and IV, (it J+ 1 — uf) is given by (5.3.11) 
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while ( Uj —Uj_i) is given by (5.3.8). The product, (u, x u, xx ), goes through zero here 
because u, xx goes through zero. The entropy production rate factors to 

A x(P t )j = (uj+ 1 - Uj - Suj-i + tty_ 2 )(— tty + i + 3 Uj - 3uj_i + Uj- 2 ) (5.3.15) 

The first factor is a first derivative formula while the second factor is a third derivative 
formula. The Taylor series expansion shows 

(P,)j = ~cu, x u, xxx Ax 3 + 0( Ax 4 ) (5.3.16) 

This shows that such transitions (inflection points) produce entropy at a lower rate than 
the basic scheme. 

Between regions II and III, u(x) reaches a maximum and u, z = 0. The computation 
of Uj + i may involve either branch of (5.3.6), depending on the relative magnitudes of the 
intervals on either side of the peak. 


If the interval to the left of the peak is larger in magnitude than the interval to the 
right of the peak, the transition point will be the first point to the right of the peak. 
At the transition point (j), is computed by the usual upwind formula (5.3.5) and 

(ttj+i — Uj) is given by (5.3.7). On the other hand, «j_i is computed by the symmetric 
formula 

Uj_x = -(uj + uj-!) (5.3.17) 

and (uj — Uj_j) is given by (5.3.12). Notice that (iiy + i — Uj) — ( Uj — Uy_i). Conse- 
quently, no entropy is produced by such transitions. 


If the interval to the left of the peak is smaller in magnitude than the interval to the 
right of the peak, the transition point will be the peak itself. At the transition point ( j ), 
Uj + i is computed by the second branch of (5.3.6). The absolute value sign is needed in 
this case and «y+i is given by 


“j+i = + 


(5.3.18) 


A simple manipulation gives 


w j) — t\( u i u i-i) 


(5.3.19) 


On the other hand, Uj_i is computed by the symmetric formula (5.3.17) and, (uy — Uj _. |) 
is given by (5.3.12). Notice that (uy+i — uy) = — (uy — Uj_±). Since these differences are 
squared in (5.3.1), no entropy is produced by these transitions either. 


To summarize, the scheme defined by (5.3.6) is second-order accurate from the point 
of view of entropy production rate errors. There is no entropy production at extrema (a 
property not shared by the standard MINMOD limiter). More importantly, the entropy 
production rate for this scheme is nonnegative on a cell-by-cell basis. 
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5.4 An Explicit Scheme Which Satisfies a Cell Entropy Inequality 

If the first-order upwind differencing scheme of Section 5.2 is advanced in time using 
explicit Euler time advance, the result satisfies the fully discrete entropy inequality for 
the wave equation provided v = < 1. This well-known result can be shown by direct 

substitution. The semi-discrete scheme is just 

(«,t)j + -^( u i ~ «i-i) = 0 (5.4.1) 

This leads to the fully discrete scheme, which, after multiplying by A t is 

u n+l _ u n + u ( u n _ ^ ) = Q (5.4.2) 


The fully discrete entropy inequality is 

c n +i _ cn F n , , — F n , 

( p\. | 

~ At + Ax 
Substituting the definitions (5.1.6) gives 

= - [(«” +1 ) 2 - («” ) 2 ] - ''[( u 7) 2 - («?-i) 2 ] 


(5.4.3) 


(5.4.4) 


The only part of (5.4.4) that is not directly computable is (u” +1 ) 2 and this can be found 
by solving (5.4.2) for u” +1 and substituting. When this is done the expression for entropy 
production rate reduces to 


A t(P.)j = ^(1 - v)(u] - uj.,) 2 (5.4.5) 

This is a positive quantity only if 0 < u < 1. This corresponds to the von Neumann 
stability range of the scheme. 


It is possible to produce a more accurate scheme if a more restrictive time step range 
can be accepted. In general the production of entropy in a single step for the wave equation 
using explicit Euler time advance is given by 



(5.4.6) 


which is derived from (4.2.8) and (5.1.8). The right-hand side must be positive to satisfy 
the second law. The last term makes this difficult, especially if v is large. Equation (5.4.6) 
can be restated as 

A t(P,)j = v(u] - ) 2 [(fl + 1) ((1 - u) - R(1 + *))] (5.4.7) 
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Figure 5.2. — Only schemes in the shaded region satisfy a 
cell entropy inequality. 


where 

(5.4. 

Assuming u > 0 (the case where v < 0 is a trivial extension), it is clear by inspection of 

(5.4.7) that ( P t )j > 0 when 

-1 < R< — - (5.4.9) 

The two inequalities are satisfied within the shaded region of figure 5.2. 

Notice that for CFL numbers ( v ) greater than one, only negative values of the ratio 
can be admitted. If ti" = then consistency suggests that = u". Referring to 

(5.4.8) this implies that R = 0. Thus, it is not generally possible to construct schemes 
having only negative values of R. 

At a CFL number of 1, the (5.4.9) reduces to = Uj, which is described above. 

For CFL numbers less than or equal to one, it is possible to construct consistent schemes 
which satisfy a cell entropy inequality. Since explicit Euler destroys entropy, the second- 
order scheme (5.3.6) described in Section 5.3 must be modified to produce more entropy. 


R = -!±i - 

u" — u? , 
J j-\ 
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This can be done as follows: 


u j+\ 


= < 


u i + “ *')(«i+i ~ u i) 


u i + jC 1 ~ v ) 


Uj - 

Uj-1 

u j+l 

-Uj 


if |u;+i - uj\ < \uj - Uj-i | 

(uj+i — u j) otherwise 


(5.4.7) 


The modified scheme which results will not be second-order accurate unless At is of 
the same order as Ax 2 . In that case (1 — u) is essentially equal to one and the scheme 
reverts to one previously analyzed (5.3.6). Unfortunately, such a time step restriction is 
generally too severe. Nevertheless, (5.4.7) is much more accurate than (5.2.1) for any 
stable value of u. 


For transient problems such as this, the spatial and time discretizations are insep- 
arably coupled. Any entropy produced by either will stay until it convects through an 
outflow boundary or until the calculation stops. It can never be removed by a scheme 
which satisfies a cell entropy inequality. 

5.5 A Stable TVD Scheme Which Violates the Cell Entropy Condition 

To construct a scheme and analyze it for compliance with the second law requires 
two things, a time advance scheme and a space differencing scheme. Space differencing 
schemes can be expressed as algorithms for computing tt J+ i. Many such schemes can be 
characterized by the formula 


«i+i = u i + ^i+*(«i+i - «j) (5.5.1) 

This requires an algorithm for computing <f>j+x • Usually 4>j+i is constrained to fall 
between 0 and 2, which is equivalent to constraining i to lie between uj and Uj + *. The 
first-order upwind scheme described above can be characterized by the formula <j>j + x = 0, 
while central differencing is characterized (for the wave equation) by 4>j+i = 1. Somewhat 
more subtle is the description of a second-order upwind scheme for which 



l 

u j+i ~ u j 


(5.5.2) 


The right-hand side is often shortened to rj. This suggests that more complicated 
schemes might be represented as a plot of i(rj). Several schemes which can be char- 
acterized in this way are plotted in figure 5.3. 
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Figure 5.3. — Many common difference schemes can be de- 
fined in terms of <f>(r) . 


A simple scheme which is stable, but violates the entropy condition, is 



r 0 

< 2 Tj 

. Tj + 1 


if rj < 0 

if Tj > 0 


(5.5.3) 


The corresponding semi-discrete scheme is spatially second-order accurate. When used 
with explicit Euler time advance, (5.5.3) is not only stable (conditionally), but Total Vari- 
ation Diminishing (TVD) (ref. 25). There it was also shown analytically that the scheme 
violates entropy. The next section will show results obtained with this scheme. 


5.6 Results 

This section shows the results of computational experiments made with schemes de- 
scribed above. Two initial waveforms were used. One was the square pulse of unit height 
and a width of 10 grid points. The other was a single period of a sine wave with the 
same amplitude and width. To compare schemes with different stability limits, and to stay 
away from the perfect shift that occurs for some schemes when v = 1, all schemes were 
implemented with v = The results are shown after 50 steps. During that time the pulse 
has propagated 25 grid points. 

The first scheme shown is the entropy violating scheme of Section 5.5, defined by 
(5.5.3). Remember that this scheme is spatially second-order accurate, conservative, stable, 
and TVD. What can possibly go wrong? In figure 5.4 we see that it does an excellent job of 
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Figure 5.4. — This scheme satisfies a TVD condition, but 
does not satisfy a cell entropy inequality. Both pulses have 
traveled 2.5 pulse widths. Notice the tendency to produce a 
square wave from smooth initial data. 


convecting a square wave. There is no noticeable amplitude loss, and the shape is virtually 
unchanged. The area under the pulse is conserved. 

The right half of figure 5.4 shows why the square wave is captured so accurately. Many 
initial pulses will eventually become square waves. This behavior reflects the decreasing 
entropy (increasing L 2 norm) within the domain. Because of the TVD condition, the L\ 
norm is bounded. This apparent contradiction is resolved only when u(x) is piecewise 
constant as it is for the square pulse. Such solutions have zero entropy production. The 
square pulse minimizes the total entropy in the domain, subject to the conservation law 
and the TVD constraint. 

The second scheme is the simple first-order upwind scheme, (5.2.1). This scheme is 
first-order accurate in space and time, conservative, and satisfies a cell entropy inequality. 
In figure 5.5 we see a substantial loss of amplitude and smearing of sharp corners that this 
scheme is famous for. 

The third scheme shown (5.4.7) also satisfies a cell entropy inequality. It is slightly 
less dissipative than the first-order upwind scheme, but remains first-order accurate in 
time due to explicit Euler time advance. This scheme is very similar to the limited Lax- 
Wendroff schemes seen in early TVD papers (ref. 11). Only the limiter is changed slightly 
as previously noted. In figure 5.6 we see that this scheme is a big improvement over the 
previous one, at least for this equation. This improvement reflects the fact that the entropy 
produced by the spatial differencing approximately cancels the entropy destroyed by the 
temporal differencing. This strategy is much less effective in systems that do not have a 
single, constant, wave speed. 
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Figure 5.5. — These results are from a standard first-order 
upwind scheme which uses explicit Euler time advance. This 
scheme satisfies a cell entropy inequality. 




Figure 5.0. — This scheme is significantly more accurate than 
the standard first-order upwind scheme. It uses explicit Elder 
time advance and satisfies a cell entropy inequality. The accu- 
racy improvement comes from a partial cancellation between 
space differencing and time differencing errors. 


A scheme which is second-order accurate in space (5.3.6) was discussed in Section 5.3. 
This scheme satisfies a semi-discrete entropy inequality. Advancing in time using implicit 
Euler time advance preserves this inequality. Since this is a linear problem and since the 
limiter is piecewise differentiable, there is little difficulty in implementing this scheme. The 
results for this scheme are shown in figure 5.7. 

The results are disappointing. The errors from the time advance dominate the calcu- 
lation. The total entropy produced is not too different from that produced by a first-order 
upwind scheme. For problems with a steady state, the time advance errors would eventu- 
ally disappear, leaving a second-order accurate solution at steady state. 
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Figure 5.7. — This scheme is second-order accurate in space. 

First-order accurate, implicit Euler time advance is used. The 
time advance errors dominate. 

To improve the time accuracy, a modified Crank-Nicolson time advance may be used. 
This was described in Chapter 4. Such a scheme eliminates the temporal dissipation in 
this case. The improvement can clearly be seen in figure 5.8. 



Figure 5.8. — This scheme uses a modified Crank-Nicolson 
time advance to achieve second-order time accuracy. Spatial 
differencing is the same as for figure 5.7. 


A single number describing the accuracy of all the schemes is total entropy produced 
during the time of integration. This is easily computed, it is just the sum of u 2 over all 
grid points for the initial condition less the equivalent sum after 50 steps. For this simple 
problem, the correct answer is exactly zero. A table summarizing these results is given 
below 
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Table 5.1. — Entropy Production For Various Schemes 


Figure 

Square Pulse 

Sine Pulse 

5.4 

0.48 

-1.87 

5.5 

3.90 

1.98 

5.6 

2.21 

0.44 

5.7 

5.40 

2.26 

5.8 

2.83 

0.86 


The clear winner is the first-order scheme of figure 5.6. Such a scheme is unsuited for 
problems having a steady state, since the time step affects the computed solution to such 
problems. Furthermore, the time step is limited to a CFL number of 1 in any event. 

The Crank-Nicolson scheme of figure 5.8 has no such problems. It is formally second- 
order accurate in space and time. The spatial accuracy is independent of the time step, 
which is limited only by the nonlinearity of the problem (see Section 4.5). Most impor- 
tantly, it satisfies a cell entropy inequality. In the next chapter, this scheme will be adapted 
to nonlinear problems. 


54 











Chapter 0. NONLINEAR SCALAR EQUATIONS 


6.1 A Suitable Entropy for Burgers’ Equation 

This chapter is concerned with the solution of nonlinear scalar equations of the form 

u, t + f„= 0 (6.1.1) 

where / is any nonlinear, continuous function of u. In particular, Burgers’ equation 

= 0 ( 6 . 1 . 2 ) 

will be explored. As previously shown, a suitable entropy function for any scalar equation 
is 

S = -u 2 (6.1.3) 

This satisfies the convexity condition. 

Syuu = — 2 < 0 (6.1.4) 

For Burgers’ equation, / = |u 2 . The compat ability condition 

F, u = S, u f, u (6.1.5) 

determines the corresponding entropy flux function F. This is 

F = -|u 3 (6.1.6) 

as may be verified by substitution. Thus, a suitable entropy pair for Burgers’ equation is 
given by equations (6.1.3) and (6.1.6). 

6.2 A First-Order Scheme 

Construction of first-order schemes which satisfy the second law was discussed in 
Section 3.4. In this section such a scheme is constructed for Burgers’ equation. 

In Chapter 4 it was shown that satisfaction of the semi-discrete entropy inequality 
implies satisfaction of the fully discrete inequality when certain time advance schemes are 
used. The semi-discrete form of (6.1.1) is 

At 

U,t+ A t 6 - 2 - 1 ) 
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where /i+i = /(«j + i ). The choice of ti J+ ^ is subject to the constraint that the second 
law is satisfied in each cell. An expression for this constraint was derived in Section 3.4: 

A = [-(S,.)i(/ j+J - fj) + (F 1+i - F,)] 

+ [-(«,.);(/; - fj- } ) + (Fj ~ Fj. j )] > 0 (6.2.2) 

where = F(uj + ± ). The terms in brackets define the right and left half cell entropy 

production rates, Axj(P+ )j and Axj(P~)j respectively. Using these definitions (6.2.2) 
can be written as 

A*j[(^)j+(A')<] >0 (6.2.3) 

The algorithm proposed in Section 3.4 is to choose tt J+ ± in stich a way that (P+ )j 
and (P“)j are individually positive. This requires that 

(P+h > 0 and (Pf) i+1 > 0 (6.2.4) 


It is not immediately obvious that a value of 1 exists which satisfies this requirement. 
The choices of v-j + i may be parameterized in general by the formula 


"i+i = u i + j “ u i) 


(6.2.5) 


where the parameter takes on values between 0 and 2 as takes on values 

between Uj and «j+i. If 0 < <f>j+i <2, and Uj+i is chosen according to equation (6.2.5) 
the following inequality holds. 


(«j+i - «j+i )(«j+i - «i) > 0 (6.2.6) 

This inequality is equivalent to (3.4.12) and allows the use of inequality (3.4.8) and equation 
(3.4.7). 

Sometimes a single example is more illuminating than a page of algebra. In figure 6.1, 
the values of Axj(P+ )j and Axj +1 (P~)j+i are plotted as functions of <f>j+i for a strong 
compression, uj = 1 and Uj + 1 = |. 

When 4>j+i = 0, = Uj and (P+ )j = 0 by inspection of equation (6.2.2). Simi- 
larly, when 4>j + i = 2, 1 = Uj + j and (Pr)j+i = 0- Figure 6.1 confirms this behavior 

for a particular case. 

Also in figure 6.1, it is clear that, for any value of <f>j+x in the range plotted, the 
slopes of the two curves have the same sign. This is a trivial extension of inequality 
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Figure 6.1. — Effect of the interpolation parameter <f> on 

semi-discrete cell entropy production rates for Burgers’ equa- 
tion during a strong compression. 


(3.4.12). Inequality (3.4.8) shows that the product of the two slopes will be positive if and 
only if (6.2.6) holds. This feature will be exploited more heavily in the next section. 

In addition to having the same sign at a given value of <f>j+i the two curves have no 
extrema on the interval if Uj and uj + 1 have the same sign. Extending (3.4.6) through the 
chain rule and simplifying gives the formula 

^ = — — (6.2.7) 

j 

which is valid for only Burgers’ equation. The identity S, uu = —2 was used in the simplifi- 
cation. Equation (6.2.7) shows that, throughout the range of <j>j + 1 shown, the slope of both 
curves is opposite in sign to the value of u J+ x. In this case u J+ i is positive throughout 
the interval so both curves have negative slopes everywhere. 

As long as Uj and Uj+i have the same sign, the choice 

<t> j+ i = 1 - sign(uj) (6.2.8) 

will always satisfy the inequalities of (6.2.4). These, in turn guarantee a cell entropy 
inequality. 

If uj and Uj+i have opposite signs, the two curves will contain extrema. In figure 6.2, 
for example, Uj — — 1 and Uj + 1 = 1. 

Notice that, despite the extrema, the two curves have slopes of the same sign for any 
given value of <f>j+i- This implies that the extrema will be at the same point on the two 
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Figure 6.2. — During a sonic expansion, both cells experience 
a maximum entropy production rate when = 0. In this 

case that occurs when 4>j+i = 1. 

curves. This is a case where neither of the endpoints satisfy (6.2.4). One point which will 
is the extremum; — 1 in this case. For sonic points such as this, the extremum occurs 

when 

=0 (6.2.9) 

It is only at sonic points such as this that the first-order scheme given here differs from the 
usual first-order upwind scheme. Incidentally, these are called sonic points by analogy to 
the one-dimensional Euler equations. They are identified by a change in sign of the wave 
speed in a region where the characteristics diverge. 


For Burgers’ equation, (6.2.9) is equivalent to = 0. Often, it is easier to imple- 
ment the formula 

4 i+ 1= (6.2.10) 

3 2 Uj+l-Uj 

because it specifies <f>j+i instead of . Substituting (6.2.10) into (6.2.5) yields ttj+A = 0. 


In figure 6.2, which depicts a sonic point, the two curves contain maxima. In the 
case of a shock, the extrema are minima. Figure 6.3 depicts such a case, with Uj = 1 and 
Uj + 1 = —0.45. 

The minimum for both curves is at 4>j+i ~ 1.4. Because they each have a single 
extremum, both curves will reach their maximum values at one of the endpoints. One 
choice which always satisfies (6.2.4) is 

= 1 - si « n ( ?, ’ + 2— ' ) ( 6 ' 2 - U ) 
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4>y* 

Figure 6.3. — A rapid compression through a moving shock. 


This evaluates to one of the endpoints, 4>j+% — 0 for this example. In cases where the 
shock is nearly stationary, roundoff error may cause the wrong endpoint to be selected. 
Fortunately, in such cases, either endpoint will satisfy (6.2.4), so this is not a problem. 


by 


In summary, a first-order scheme which satisfies a cell entropy inequality is defined 




1 - sign(uj) 


sonic points ( Uj < 0 < Uj+ 1 ) 

shocks (uj > 0 > Uj+i) 
elsewhere (ujUj+i >0) 


( 6 . 2 . 12 ) 


The resulting scheme is exactly equivalent to a standard first-order upwind scheme 
except at shocks and sonic points. It is precisely at these points that entropy production is 
important. In the case of sonic points, entropy must not be destroyed in a numerical expan- 
sion shock. In the case of shocks, entropy must be produced as required by physics. Any 
numerical scheme without these properties will generate qualitatively incorrect solutions. 

6.3 A Second-Order Scheme 


The scheme devised in Section 6.2 is consistent and stable, but only first order ac- 
curate. Substantially more entropy is produced than is required by the physics. This 
section explores a variation which is second-order accurate. The technique was outlined in 
Chapter 3 and demonstrated in Chapter 6 on the linear wave equation. 

A scheme which is consistent and second-order accurate is 

= 1 («- 3 -l) 
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This is not quite the same as central differencing because f(u) is nonlinear. Although 
this scheme does not (in general) satisfy a cell entropy inequality, it may be possible to 
construct a scheme which does, in such a way as to approach (6.3.1) in the limit as the 
grid is refined. Such a scheme would be second order accurate in this limit. 

The cell entropy inequality for the j th cell is given by 

(p.)i = ( P. + )i + (>.-); > 0 (6.3.2) 

Where the first-order scheme of Section 6.2 required that these terms be individually 
positive, the second-order scheme requires only that their sum be positive in each cell. 

The entropy production rate for the scheme defined by (6.3.1) might be written 

(P.)i = W)i + WTi (6-3.3) 

where (Pf~ )j is the value of (P+ )j which corresponds to = \(uj+i +uj). The entropy 
production rate (P 3 )j might be positive or negative, depending on the data. 

Now consider schemes for which 

Wi>WTi and (P;)i>(PT)i (6.3.4) 

Under this convention the inequality 

Wi + > ( P*)i + WTi (6.3.5) 

holds. If the right-hand side is positive, the left-hand side is also positive and the entropy 
inequality (6.3.2) is satisfied. This is helpful since the right-hand side of (6.3.5) has only 
a single free parameter, <t>j+i, which determines u J+ ± (and ultimately (P 4 + )j). Unfortu- 
nately, satisfaction of (6.3.5) may not be possible; for example, if tty = Uj+i and (Pr)j < 0. 
This necessitates modifying (6.3.4) as follows: 

Wi > min (-WT,0) . (6.3.6) 

(P.-)i > inn (-(P, + )/,0) (6.3.7) 

In the last section it was shown that can always be chosen in such a way as to 

make (P?)j and (P~)j+i nonnegative. This is sufficient to show that relations (6.3.6) and 
(6.3.7) can always be satisfied. In combination with (6.3.4), these are sufficient to satisfy 
a cell entropy inequality. 
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Since both (6.3.6) and (6.3.7) must be satisfied for each cell, the choice of 4>j+i must 
satisfy am additional inequality which is just (6.3.7) shifted by one index 

(/>.-),+! > min (-(/>+),'+., o) (6.3.8) 

This is illustrated graphically in figure 6.4. As in figure 6.1, Uj = 1, and Uj+ \ = 
0.5. Since (P~)j+i > 0, condition (6.3.8) will automatically be satisfied. The horizontal 

chaindashed line indicates the value of — (P,”)i- This value depends on Uj and Uj-i 
(uj-i =1.25 in this case) and is solution dependent. From the figure, one can see that a 
vadue of 4>j+i » 0.75 would be chosen. This satisfies the inequalities (6.3.4), (6.3.6), and 

(6.3.8) amd thus meets the cell entropy inequality for the j th cell. If the vailue of — (P*“)j 
were lower (perhaps -0.1), a value of 4>j+i = 1 would suffice. This is a much less dissipative 
method than the first-order scheme of Section 6.2. That scheme would have chosen the 
vadue <f>j+i = 0. 



Figure 6.4. — Lowering the vadue of <f> to 0.75 guarantees a 
positive value of (P s )j- The value of (P a )j+i is also increased 
slightly, a side effect of the technique. 


6.4 Practical Considerations 

The scheme described in Section 6.3 requires computing the intersection point between 
a horizontal line and a curve. This was done graphically in figure 6.4; in practice a modified 
Newton’s method can be used. There is also an elegant approximation that works well for 
Burgers’ equation. 
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The Newton’s method approach uses the expression for A Xj(P+)j which is, from 

( 6 . 2 . 2 ) 

A*i(P. + )y = 1 - fi) + (F i+i - Fi) (6-4-1) 

Using the expressions given in Section 6.1 for this becomes, for Burgers’ equation 

^2 U,;+ * 2 ~ 2 U 3 U i + J 3 3 Uj3 ) (6.4.2) 

This expression factors 

Ax j(P> + )i = ~ i - + |(«;+± ~ “i)) ( 6 - 4 - 3 ) 


Equation (6.2.5) expresses u J+ i in terms of Uj+i and <f>j + 1 . This allows computation of 
the gradient which is required for Newton’s method. 


Ax j - u if < i > i+\ u i+\ (6.4.4) 

This expression for the gradient has three zeroes which must be guarded if Newton’s 
method is to be reliable. In the first case, if «j+i — Uj = 0, the choice of 4>j+i is arbitrary 
and has no effect on the choice of Uj+i- In the second case, 4>j+i — 0, which can occur 
during strong compressions or shocks. The third case, u J+ i = 0, occurs at sonic points 
and shocks. The value of 4>j+i at which this occurs is 



2uj 

ttj+1 U j 


(6.4.5) 


In the first case, the choice <f>j + 1 = 1 is made. In the second and third cases, the 
change suggested by one step of Newton’s method is limited to 90% of the distance to the 
nearest extrema. Finally, the magnitude of the change is limited to keep 4>j+i between 0 
and 2, in case the initial guess is very close to an extremum. 


In practice, this works quite well and converges in about four iterations to six decimal 
places. It works equally well for (P~ )j for which the expression is 


Ax j(P. )j = (uj - - ^( Uj - «-,_i)^ 

and for which the gradient is 


Axj ^izr = - u ;-i ) 2 ( 2 - 


(6.4.6) 


(6.4.7) 
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At sonic points, it is possible that both (jP , + )j and (P, )j are negative. This can result 
in the null stencil, /i+i = fi-\ = ft = which implies ( u, t )j = 0. This only occurs if 
the sonic point happens to lie on a grid point, in which case is supposed to be zero. 

At shocks (P+)j-_i and (P a “)i ma y both be negative. This results in two different 
solutions for <t>j~± • The value farthest from 1 is selected naturally by Newton’s method. 

An approximation is often made to be computationally efficient. Referring to equation 
(6.4.3), if |(u J+ 1 — Uj) is small compared to uj, it can be neglected. In that case (6.4.3) 
reduces to 

A *i(A + )i » -«i) 2 ( 6 * 4 * 8 ) 

Equation (6.4.8) can be solved with only a square root, but even this can be avoided. If 
(P“)j is approximated in the same way as (P+ )j , the result is 

Az i(-P» )i ~ — ~ u i-0) 2 (6.4.9) 

Summing (6.4.8) and (6.4.9) gives an approximation for the entropy production rate. It is 
this approximate entropy production rate which can efficiently be made positive. Evalu- 
ating (6.4.9) at = 1 gives an approximation for (P,~)j. 

Axj (P,~ )j « ij(uj - Uj-i ) 2 (6.4.10) 

Using these approximations, it is easy to satisfy the sufficient condition (6.3.5). With very 
little algebra it can be expressed (for uj > 0) as 



(6.4.11) 


The right-hand side is just the easily computed r ; 2 . An equivalent expression is 


^ l r il ( 6 - 4 ‘ 12 ) 

The convention (6.3.3) requires that <f>j + i < 1 when uj and Uj+ 1 are both positive. So in 
this case, a sufficient constraint for satisfying the approximate entropy inequality is 


< min(l, | r j | ) 


(6.4.13) 


It is worth noting that this limiter is identical to the limiter used for the wave equation. 
In fact the approximation used amounts to assuming a constant wave speed between « •_ i 
and This is essentially the same approximation (but not quite the same limiter) 

that is used in TVD schemes. 
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Expression (6.4.10) is only valid if uj and «j+i are positive. If both are negative, the 
equations are different but the steps are the same. The limiter in this case turns out to 
be, 


+i+\ ± )) 


(6.4.14) 


This expression can also be derived by reversing the coordinate system (so that uj and 
Uj+i are positive) and then using (6.4.12) directly. 


The approximation that Uj+ 1 — Uj is small compared to uj is especially bad at shocks 
where the former is large and the latter is small. For this reason the first-order scheme is 
used at shocks. 


6.5 Results 

If the approximation of a locally constant wave speed is made, this scheme becomes 
formally TVD. In the numerical experiments which follow, the locally constant wave speed 
approximation (6.4.13) is used. To check the approximation, <f>j+i is also computed with- 
out the approximation by using Newton’s method. The two values of are compared. 

Where appropriate, the effect of any discrepancy on the solution is noted. 

The initial condition chosen here is a two-period sine wave of unit amplitude. This 
problem contains sonic points and develops shocks. Given that the domain is of length 
47 t, the shocks will first occur at time t = 1, rapidly strengthening until time t = j . At 
the shock, the exact solution is double valued. After t = j, the shocks diminish steadily 
in amplitude while the solution maintains its double sawtooth shape. The steady state 
solution is u = 0. 

In figure 6.5 the initial condition is plotted. The open diamonds represent the com- 
puted solution, while the solid line represents the exact solution. The boundary conditions 
are Dirichlet, but periodic boundary conditions would give the same results in this case. 

Figure 6.6 shows the initial plots of <f>(x). In this case, the solid line represents the 
values computed with Newton’s method (labeled as exact) and the diamonds represent the 
values computed with (6.4.13) and (6.4.14) (labeled as approximate). 

Notice that in regions where (u 2 ), x > 0, the approximate limiter is just <f>j + i = 1 
which approximates centred differences. In compression regions, some dissipation is added, 
proportionately more as the compression becomes stronger. Here the method approximates 
the Warming-Beam second-order upwind scheme. In a sense, the scheme simply switches 
between a second-order upwind scheme and a central difference scheme. 

When \uj\ is much larger than |(iij +1 — «j)|, (6.4.8) and (6.4.9) are good approxima- 
tions for (Pf)j and (P“)j. In these regions, (6.4.13) ( u > 0) and (6.4.14) ( u < 0) yield 
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Figure 6.5. — Burgers’ equation will sharpen this initial con- 
dition into a two-period sawtooth wave as time passes. 



Figure 6.6. — The value <f> = 1 corresponds to a symmetric 
interpolation. Other values indicate various degrees of upwind 
bias. Two different ways of computing <j> are compared here 
for the initial condition. 

nearly the same value of <f>j+r as the much more costly Newton’s method. 

This is not true however, when Uj is near zero. Referring to figure 6.6, there is 
significant disagreement between the two methods in small regions around sonic points 
and “shocks.” These regions become smaller as the mesh is refined. 

At x = 7r and x = 3ir where u changes sign, the exact <f> is continuous while the 
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approximate <f> is obviously discontinuous. At x = 2 tt, the sonic point, the reverse occurs. 
The approximate <f> is continuous while the exact <f> is discontinuous. 


The difference between the exact and approximate values for <f> can be explained 
entirely by the nonlinear variation of wave speed across each cell. Consider a case in which 
u(x ) is linear and positive. The exact entropy condition (6.2.2) reduces to 




.1 

/ \2 

2, .1 

1 

**■* 

3 

~(u j+ x - Uj) 

Uj + -Uj) 


(6.5.1) 


If (ttj+i — Uj) = (uj — Uj_x) (i.e., no dissipation), then satisfaction of the inequality 
depends on the signs of (uj — u,-_i) and (tt J+ i —«_,). Around the sonic point, both are 
positive, requiring some adjustment of <f> to make the magnitudes different. On the other 
hand, when both are negative, the effect is beneficial, allowing the dissipation to be turned 
off altogether. 

The constant wave speed approximation obscures this effect. With the approximation, 
the entropy condition reduces to 



a) 2 -(“;+$ 


-Uj) 2 > 0 


(6.5.2) 


Satisfaction of this inequality depends only on the magnitudes of the two terms and not 
on their signs. Around a sonic point, both terms have approximately the same magnitude 
so <j) is adjusted very little, if at all. Around a shock, the first-order scheme is used. 

In figure 6.7, the semi-discrete cell entropy production for each cell is given. The sym- 
bols show the values for the approximate <f> while the solid line shows the values obtained 
from the exact <f>. These values are independent of the time step. There is very little 
entropy production initially because there is no shock. What little there is, is numerical 
and should be eliminated if possible. 

The solid line has no negative values, as planned. This shows that a cell entropy 
inequality can be satisfied in principle. The symbols, on the other hand show negative 
values near the sonic point, and near the boundaries (which are essentially sonic points). 
They also show a relatively large amount of entropy production in the cells at z — it and 
x = 27t where the first-order method is applied. The accuracy penalty exacted by the 
approximation is apparent throughout the entire compressive region. 

In carrying out the time integration, the scheme 


(6.5.2) 
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Figure 6.7. — Initial semi-discrete cell entropy production 
rates. Notice the negative values near the sonic points. 


is ideal for scalars. As explained in Section 4.4, the fully discrete entropy production of this 
scheme is the same as the semi-discrete entropy production, independent of time step. The 
scheme is very difficult to implement exactly because of linearization errors. As explained 
in Section 4.5, the following variant is much more practical 


u 


n+| _ 


1 At 


u 


n+l 




(6.5.3) 


The fluxes f*. x and /*_ , are computed with some suitable linearization so as to 

requires derivatives 


approximate / 

a*. 


n+* 

j+i 


and f n +i . Computing derivatives, such as ■ ,+ ^ 
* 2 


f -J~t - - 9 h 

such as a 2 . This is where equations (6.4.12) and (6.4.14) exhibit a tremendous compu- 
tational advantage. The approximate definition of </>j+i is piecewise linear and derivatives 
can be obtained analytically. By contrast, if <f>j + 1 is computed with Newton’s method, 
numerical differentiation is needed. For this reason, the approximate <f> values are used to 
advance in time. The maximum CFL number used will be 0.5. 


In figure 6.8, the solution is plotted at time t « 1.005. The sharpening of the solution 
into a shock is quite apparent. There are no oscillations anywhere and the computed 
solution is tracking the exact solution quite well. There is a slight loss in amplitude at the 
peak values, perhaps due to the excess dissipation in the compressive regions. 

In figure 6.9, the values of (f> are shown at the time of figure 6.8. Clearly the compres- 
sive region has shrunk and become more intense, while the expansion region now occupies 
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Figure 6.8. — Onset of the shock. 


most of the domain. Notice how the solid line has begun to respond (symmetrically) to 
the shock. Furthermore, the approximate <f> is becoming more symmetric as it responds to 
what it perceives as a very strong compression next to a shock. In fact, there is no way to 
distinguish between a strong compression and a shock on a finite mesh. 



Figure 6.9. — Dissipation at shock onset. 


In figure 6.10, the semi-discrete entropy production is shown. Notice the larger values 
of entropy production that have begun to appear at the shock (the scale has changed by 
three orders of magnitude). Also notice the way in which the dissipation has followed the 
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Figure 6.10. — With the appearance of the shock comes 

significant entropy production within the shock cells. For ex- 
ample, compare scales with figure 6.7. 


solution so that essentially all of the dissipation is at the shock as required by the physics. 
Finally, the entropy production rate is higher than it should be. 

At t —■ j the peak of the sine wave has traveled exactly 12.5 grid points and the peak 
shock strength occurs. The computed solution is shown in figure 6.11. The computed 
solution is very smooth and tracks the exact solution except at the shocks. There is again 
a slight loss of amplitude. 



Figure 6.11. — This is as strong as the shock gets. After this 
it slowly decays to zero. 
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In figure 6.12, the compressive regions have shrunk to a two-cell interval in the neigh- 
borhood of the shock. The remainder of the domain contains expansion regions. Notice 
that <f>j + 1 = 1 except at the shock. This corresponds to the analytic solution in that all 
of the dissipation occurs at the shock. 



Figure 6.12. — Essentially all of the dissipation introduced 
by the numerical method is at the shock. This corresponds to 
the physical situation. 


The peak shock strength can be computed analytically. As shown in figure 6.13, the 
peak entropy production rate for the cell containing the shock is exactly | . The correspond- 
ing numerical shock strength (represented by the diamonds) is about 1.21, a significant 
error. This is due to numerical dissipation inherent in the method. Entropy production 
rates away from the shock are less than 10~ fl in magnitude, though, as previously noted, 
some of them are negative. 

Interestingly, the computed shock strength does not decay quite as fast as the exact 
shock strength. Indeed, the numerical shocks eventually become stronger than the exact 
ones. Perhaps this is due to the entropy violations observed early in the calculation. These 
influences would not be felt at the shock until about t = tt which is about what is observed. 
Note that the entropy production errors at sonic points do not jeopardize the stability of 
the computation. Since the totcil entropy produced in the domain is positive, the stability 
proof of Chapter 2 applies and the solution remains bounded. 

It is now reasonable to review the costs and benefits of the locally constant wave 
speed approximation. This approximation results in a scheme that is more dissipative 
than it needs to be during transient compressions, though the effect is not severe. It also 
results in minor violations of the second law in the neighborhood of sonic points. These 
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Figure 6.13. — The analytic solution reaches a peak entropy 
production rate of 1.33 . The shortfall indicates accumulated 
error. 


errors can contaminate the solution elsewhere but the effect is small. On the other hand, 
the computational convenience of such an approximation is very significant. The expense 
of using Newton’s method to compute <j>j+i is quadrupled by the requirement of using 
numerical linearizations. In sum, the practical advantages of the approximation more 
than make up for its theoretical shortcomings. 

The time advance scheme used provides at least an a posteriori evaluation of the fully 
discrete entropy production. In this way the importance of the linearization errors can 
be glimpsed. In figures 6.14 and 6.15, the solid line represents the fully discrete entropy 
production, that is the total entropy produced in a given cell over the time step divided 
by the size of the time step. The symbols represent the semi-discrete entropy production 
halfway through the first step, figure 6.14 corresponds to a CFL number of 2. Notice 
that the fully discrete inequality corresponds very well with the semi-discrete scheme. 
This shows that the linearization errors are not too important at this time increment. At 
smaller time steps the situation improves further. At a CFL number of | there is no 
discernable difference. 

On the other hand, figure 6.15 corresponds to a CFL number of 5 which is distinctly 
unstable. The fully discrete entropy production rates are significantly different than the 
semi-discrete values. Some violations of the second law are apparent even in the first 
step. This indicates large errors due to linearization in time. If (6.5.2) could somehow be 
implemented exactly, instead of in a linearized sense, the scheme would be unconditionally 
stable. 
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Entropy Production Rate • Entropy Production Rate 



Figure 6.14. — At CFL = 2, the semi-discrete entropy pro- 
duction rates are in good agreement with the fully discrete 
rates. This indicates that the time linearization errors are 
small. 



Figure 6.15. — At CFL = 5, the semi-discrete entropy pro- 
duction rates are no longer in agreement with the fully dis- 
crete rates. Negative entropy production rates are evident in 
several cells. This indicates that the time linearization errors 
are significant. To reduce them requires a smaller time step or 
iteration within a time step. 


Burgers’ equation provides a good forum for testing approximations and implementing 
schemes. In this chapter, the problems of shocks and sonic points have been explored and 
the approximation of constant wave speed has been investigated. The effect of a finite time 
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step in the presence of linearization errors has also been investigated. Most importantly, it 
has been shown that a semi-discrete cell entropy inequality can be systematically satisfied, 
at least for the scalar case. 

In the next chapter these ideas will be extended to systems of equations. This will in- 
troduce several additional complications, among them: multiple wave speeds, definition of 
<j> and linearization difficulties. The fundamental principles remain unchanged. Satisfaction 
of a fully discrete cell entropy inequality (plus consistency and conservation) guarantees a 
stable and accurate solution. No other paradigm can promise this. 
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Chapter 7. ONE-DIMENSIONAL GASDYNAMICS 


7.1 The One-Dimensional Euler Equations 

This chapter is concerned with the numerical solution of the one-dimensional Euler 
equations. These were introduced in Chapter 2. For convenience, selected portions will be 
repeated here. 

In one dimension the Euler equations can be written 


where 


9>t + /,* — 0 


p ' 


pu 

pu 

and / = 

pu 2 + p 

e 


(e + p)u 


The quantities p and it represent the local mass density and the local velocity, respectively. 
The quantity e represents the total energy density. An expression for the pressure p is 


p=( 7 -l)(e--/>« 2 ) 


(7.1.3) 


A suitable entropy pair for inviscid, one-dimensional gasdynamics is 


where 


S = ps 

(7-1.4) 

F = pus 

(7.1.5) 

II 

cw 

(7.1.6) 


As shown in Section 2.3, this pair agrees with the physical viewpoint and satisfies the 
convexity and compatability conditions. When the problem is made more complicated 
by the addition of viscosity, heat conduction, and three-dimensional domains, the entropy 
does not change at all and the entropy flux changes only slightly. Although this chapter 
(and indeed, this entire work) does not proceed beyond the simple case of one-dimensional 
inviscid flow, there is little conceptual difficulty in doing so. This is discussed further in 
Chapter 10. 
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7.2 Entropy Variables 


The main difference between equations (7.1.1) and equation (6.1.1) is that q is a vector 
and u is a scalar. This introduces a number of ambiguities which need to be discussed. 
The first issue to be decided is the choice of independent variables. A common choice is 
the conservative variables. 


9 = 


P 

pu 

e 


(7.2.1) 


These are called the conservative variables because domain integrals of them are 
conserved physically (except for boundary conditions). They facilitate construction of 
schemes which conserve mass momentum and energy numerically. 

Other variables are equally valid. For example, primitive variables are often used 
because of the ease with which they can be measured experimentally. These are 

(7.2.2) 



There is a one-to-one correspondence between the conservative variables and the prim- 
itive variables. That is, each set can be computed from the other. If computational expense 
is not an issue, any scheme which can be implemented in one set of variables can be im- 
plemented in the other, using two conversions per step. 

The issue of conservation is independent of the choice of variables. A scheme conserves 
mass, momentum and energy if the computed fluxes are the conservative ones 


/ = 


pu 

pu 2 + p 

(e + p)u 


(7.2.3) 


and if these are evaluated on the control volume boundaries. These fluxes can be evaluated 
with equal ease using either the conservative or the primitive variables. 


From this point of view, any set of variables will do, provided there is a one-to- 
one correspondence with the variables of (7.2.1). This is sufficient to ensure that the 
conservative fluxes can be computed and conservative schemes constructed. The issues of 
consistency and satisfaction of the second law remain. 


In Chapter 2 the definition of consistency was given. If information at several grid 
points is used to compute a flux, and if the values at those grid points are all the same, then 
the computed flux must match the flux at the grid points. This is a fairly easy criterion 
to meet. In this work, it is met by requiring that qj+i be between qj and qj+ 
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In the case of scalars, it is clear what “between” means. The value u J+ i is said to be 
between Uj and Uj+i if 

(u i+ , - ti i+ i)(u i+ i - ttj) > 0 (7.2.4) 

In the case of systems, the definition is less obvious. In this chapter, fy+i is said to be 
between qj and qj + \ if 


V dq >+l ) V ) 


(7.2.5) 


As shown in Section 3.4, (7.2.5) reduces to (7.2.4) in the case of scalars. 

The expressions for (P+)j and (P~)j+i are (as derived in Section 3.4) 
A *i(P. + ), = - fj) + (F j+i - F,) 

A*i+i(-P.”);+i = -(-^.?)j+l(/j+i - /j+i) + (Fj+1 - Fj + ±) 

and the gradients used in (7.2.5) are 


A*j+I = [( S -«W " < S «)i+s] (/*)j+ 


(7.2.6) 

(7.2.7) 

(7.2.8) 

(7.2.9) 


Suppose that qj+\ = qj. By inspection, the two gradients have the same magnitude 
but opposite sign. Condition (7.2.5) only holds if both gradients vanish, that is if (S, q )j = 
= (S>q)j+ 1* 

Remarkably, this is sufficient to imply that qj = 1 = qj+i. Let S, q define a new 

set of dependent variables, v, that is, from (2.3.13) 


v — S, q T — [s — (7 + 1) + 




irzlkl 


p 


]' 


(7.2.10) 


Since v, q = S m , and by convexity S m < 0, it is clear that v, q ^ 0. In general, this is 
sufficient to ensure a one-to-one mapping between q and v. In particular, for the one- 
dimensional Euler equations, v can be computed from q using (7.2.10) and q can be com- 
puted from v by computing u,s,p,p, and e in the following sequence: 

t>2 


u = — - 


V3 

1 

s = Vi + 7 + -Ul>2 


p = 

p = 
e = 




(7-1) 

[v 3 exp(s) J 

(7 ~ ])P 

V3 

P l-i 

(P T) + 2 pu 


(7.2.11) 

(7.2.12) 

(7.2.13) 

(7.2.14) 

(7.2.15) 
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where the so called “entropy variables,” t>i, v 2 , and v 3 are the elements of v as shown in 
(7.2.10). This transformation is due to Mock (ref. 26). 

To summarize, if qj = qj+i, and condition (7.2.5) is met, then qj = = qj+i- It 

follows that fj = fj + 1 = fj+i- Thus (7.2.5) is sufficient to assure consistency in the sense 
of Lax. Careful examination of (7.2.8) and (7.2.9) (the factors of (7.2.5)) shows that they 
are easily written in terms of entropy variables (see Section 7.4). Thus, it is useful to think 
in terms of the entropy variables v when constructing consistent schemes. 

7.3 A Remarkable Identity 

In constructing schemes that satisfy an entropy inequality, the basic approach for 
systems is to diagonalize the equations locally. This can be done without approximation 
for the one-dimensional Euler equations, mostly because of the existence of a remarkable 
identity. 

Recall from Section 2.3 that the convexity condition was demonstrated for these equa- 
tions using the factorization 

s,„ = -(r , ) T P M ) (7-3.1) 

It is well known that the flux Jacobian also has a factorization 

/„ = YAY 1 (7.3.2) 


where 



u + c 


(7.3.3) 


contains the eigenvalues of f, q and the columns of Y contains the eigenvectors of /,,. It 
is well known that each eigenvector is unique only up to a multiplicative constant. It is 
possible to choose these constants in such a way that (7.3.1) holds, using the same Y as 
(7.3.2). 


The matrices Y and Y' 1 are repeated here with the correct scaling. The speed of 
sound, c = has been used to simplify the expressions. 


>> 

fit 

/3 2 {u - c) 


02 

(3 2 (u + c) 


- UC+ 

rry) 



2c a — ( 7 — l)u a 

2(7-1)u 

-2(7-1) 1 

(i) 

01 

( 7 -l)u a , uc 

Pi 

-( 7 -l)u-c 

Pi 

(*r-i) 

2/3j + 

( 7 — l)u a tic 

- 2/3a p 2 

Pi 

-( 7 -l)u+c 

02 

Pi 

(-r-D 

02 


01 ft / 2t ( r- 1) 


(7.3.4) 


(7.3.5) 
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The fact that fa is used twice is not surprising. The equations are expected to be invariant 
to a reversal of the coordinate system, so the u + c and u — c waves are scaled the same 
way. 


The surprising thing is that (S',,,)' 1 can be decomposed in terms of the eigenvectors 
of /,,. This is a new result. It has been verified that this is also true for the shallow water 
equations in one dimension (a two-by-two system) 1 , and of course, for scalar equations. So 
far, no theoretical explanation is offered. The identity itself was found and verified using 
MACSYMA. 

There are some interesting corollaries. It will be necessary in the next section to deal 
with the flux Jacobian with respect to v, the entropy variables. This is related to the usual 
flux Jacobian and its eigensystem as follows: 


/>« — f yv v q 

chain rule 

(7.3.6) 

= f yvSiqq 

definition of v 

(7.3.7) 

yay - 1 = -/ w (r 1 ) T (r I ) 

(7.3.1) and (7.3.2) 

(7.3.8) 

/,„ = -y\y t 

from (7.3.8) 

(7.3.9) 

Equation (7.3.9) is consistent with the well known symmetry of /,„. 



The term 5,„/„ also appears from time to time. Using (7.3.1) and (7.3.2) immedi- 
ately gives 

S,„f„=(Y-') T HY') (7.3.10) 

which is also a symmetric form. 

In summary, identity (7.3.1) provides a new tool with which to study hyperbolic 
systems of equations in the context of the second law. In the next section it allows the 
equations to be diagonalized without approximation. 

7.4 A Second-Order Scheme 

For this problem, construction of a first-order accurate scheme which satisfies a semi- 
discrete cell entropy inequality is nearly as difficult as constructing a second-order accurate 
scheme with this property. For this reason, the first-order scheme is left to the interested 
reader. 

In this section, a second-order accurate scheme is described, which also satisfies a semi- 
discrete cell entropy inequality. The method presented here is analogous to the one for 
Burgers’ equation, but with additional complications. The primary emphasis in Chapter 6 
was on suitable approximations to make the computation more practical. In this chapter, 
the emphasis is on rigorously satisfying the semi-discrete entropy inequality, without regard 
for the practical aspects. The use of implicit Euler time advance allows the semi-discrete 
entropy inequality to be extended to a fully discrete one as explained in Section 4.3. 

1 P.L. Roe, personal communication 
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Recall that the consistency condition is 


V A+! J v / “ 


(7.4.1) 


The gradients, previously given in terms of q ((7.2.8) and (7.2.9)), can be restated as 

a diP+)j r IT 


~ v i\ (/w);+± 

(7.4.2) 

1 T 

V i+1 - 1 (/w)j+A 

(7.4.3) 


»i+i 


Consistency only requires that (7.4.1) hold when qj = ?j+i. However, other aspects of 
the scheme require (7.4.1) to hold more generally. By inspection of (7.4.2) and (7.4.3), a 
special case occurs when Vj + 1 is chosen to be the arithmetic mean, v = |(uj + vy+i). This 
choice has the property that 


Ax Ax 

& x j-a - ** x 3+1—EZ 

dq j+ i Oq j+ i 

which certainly satisfies (7.4.1). Of course q = q(v) also does. 


(7.4.4) 


Given that q is defined in terms of v, it seems reasonable to define the gradients in 
terms of v as well. The algebra for this requires the compatability condition in terms of 
entropy variables. This is derived as follows: 


5, g /, g = /,, compatibility condition (7.4.5) 

SiqftqQtv = t> (7.4.6) 

S, q f, v = F, v chain rule (7-4.7) 

Recall that (Pf)j and (Pr)j+i are given by (7.2.6) and (7.2.7). Differentiating these with 
respect to v J+ 1 gives 

Ati ThT7 = + < F ” )>+* < 7 - 4 - 8 > 

Asj+. ^'A 1 = - (F,„ ) Ht (7.4.9) 

J+ 2 

Now applying the compatability condition in entropy variables (7.4.7) gives 

Axi &TT = + (S„h +i (M i+i (7.4.10) 

Axi+I ^P t 1 = - ( s w), +i (/..);+ j (7.4.11) 

J+ 2 
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Collecting terms gives 


dv 


j+i 

7 )j±± 

'i+\ 



(7.4.12) 


(7.4.13) 


Finally, using the definition of v (7.2.10), gives 

Axi i>prf = (v ' + i " < ’’ )T(/i " ) ;+i 

3-r 2 

d{P7)t m 


^i+4 


= ( v i+i -^+i) T (/,»)i + i 


Again, if Vj +1 = v, 


d(P7)i 


r Otfr)i+i 

Cj+, ~0^r 


(7.4.14) 

(7.4.15) 

(7.4.16) 


Ax j ^=Ar J 

which is analogous to (7.4.4). It is, however, important to explore the behavior of these 
gradients away from = v. Identity (7.3.9) is very helpful in this regard. 


/,„ = -YAY* 


(7.4.17) 


Recall that A and Y are the eigenvalues and eigenvectors of f, q (not M Using this 
identity, equations (7.4.14) and (7.4.15) may be rewritten 

KTT = ~ v i'> TYKYT ( 7 - 4 - 18 ) 

2 

A * i+1 ^SpT 1 = - (vi+l - i )Ty AyT (7 - 4 - 19) 

J+ J 

where both Y and A are evaluated at q(vj + x )• Bringing Y inside the parenthesis gives 

Axj = ~(Y T v j+ x - Y T Vj) T AY T (7.4.20) 

° v i+\ 1 

A*i+ ■ - Y T v Hi f\Y T (7.4.21) 

OV j+\ 


It is convenient to define characteristic variables. 



where Y T is always evaluated at the still unknown state, q(vj + i). For example, (7.4.22) 
implies vj = Y T vj T even though Y T is evaluated at With this definition (7.4.20) 

and (7.4.21) may be rewritten as 

Ax *Wt = _ 5 ' )TArT (7 ' 4 ' 23) 

Azj+1 %PpJ ±l = _( ij+1 _ v j+ i) t \Y t (7.4.24) 

Now, using the chain rule, the left-hand side can be rewritten in terms of gradients with 
respect to the characteristic variables. 




-vj) T AY t 


-(v i+i -v j+1 fAY T 


(7.4.25) 

(7.4.26) 


Definition (7.4.22) can be used to carry out the differentiation and show that 


dv i+ i 

= Y t (7.4.27) 

dv j+\ 


This leads to the expressions 


Ax ‘Wt yT = - ( *» i~^ TAYT 

J+ 2 

Ax »'^Pr YT = -«i + j) T Ai' T 

]+ j 

Both sides are right multiplied by (F T )‘ X leaving the expression 


(7.4.28) 

(7.4.29) 


* d(P+)j _ \T * 

^ %iT"” i+i “ 

= "(w “ 5 i+*) TA 

j+ 2 


(7.4.30) 

(7.4.31) 


Notice that the gradients are now defined with respect to v instead of v. The unknown 
quantity Vj + 1 is now interpolated from the endpoints using the formula (analogous to 
(6.2.5)) 

Vj+1 = Vj + ^(f> j+ i D {vj + i - Vj) (7.4.32) 
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where the notation 4>j+± D indicates a diagonal matrix whose entries are the entries of the 
three component vector These components, <f>\ ,<f >2 , and <j>z vary in value between 

zero and two. Consistency immediately follows if Vj+i is chosen using (7.4.32). Using 
(7.4.32), equations (7.4.30) and (7.4.31) can be rewritten as 


A Xj 






S(P.-)hi 

d ”j+k 


= - *j) T ( 21 - 4>j+i D )k 


(7.4.33) 

(7.4.34) 


When 4>j+i D = I, the two gradients are equal. As in the scalar version of this scheme, 
<f>j+i may be perturbed from this value to increase the entropy production rate, (P,)j or 

(P t )j+i- If a given change in an entropy production rate is required, the smallest possible 
change (in the sense of Euclidean length) in 4>j+i is desired. Such changes occur along the 
gradients of entropy production with respect to <f>j+x • These can be computed as follows: 
first note that (7.4.32) may be written 


5 i+i = + \iVi + 1 “ 


This implies that 


The chain rule implies that 


fty+±_ h- 

d<fi j+ x 2 




dP. dp, dv j+i 


d< t>j+\ dd j+ \ d< t>j+\ 

Using (7.4.33) or (7.4.34) for the first factor and (7.4.36) for the second yields 

Ax i !r ^ - 5 i) T ^i+i DA ( 5 i+i -i>i) D 


It 


= ~4 ( i5 i+ 1 “ i>j) T ( 21 ~ 4j+\°) A (®j+i “ Vj) D 


(7.4.35) 


(7.4.36) 


(7.4.37) 


(7.4.38) 


(7.4.39) 


The right side of (7.4.38) is a three component row vector, the first component of which 
is — |(v J+ i —vj ) Consistency requires 0 < <f>\ < 2. The only factor which can be 
negative is Ai, one of the local wave speeds and the first component of A. This implies 
that <f > i may have to be reduced in magnitude if Ai > 0, or increased if Ai < 0. Thus the 
concept of upwind differencing is recovered from the second law, for systems of equations. 
No approximations were needed for this result. Note that this result does not constitute a 
scheme; although the gradients are exact, they only give a direction and an approximate 
magnitude for perturbations to a given vector 
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The matrix <j>j + 1 D must have the same sign as 21 — <f>j + 1 D , term by term, since the 
entries are constrained to lie between 0 and 2. Since A is evaluated at fy+i = q(v J+ i ), 
both gradients must both have the same sign, term by term. 

All this information on how the entropy production rate depends on 4>j+ii is needed 
to construct a second-order accurate, consistent scheme which satisfies a cell entropy in- 
equality. A semi-discrete expression for the entropy production rate is 

A x i(P.)i = A*,( P+)j + A Xj(Pr)j (7.4.40) 


The initial guess for fy+i will be called q. It is constructed by computing i>, the arith- 
metic mean of Vj and Vj+ 1 , and then reconstructing q(t>) using (7.2.11) through (7.2.15). 
This corresponds to <f>j + 1 = 1. Using Qj+i — q, it is possible to construct (Pp)j and 


(P» )j by using (7.2.6) and (7.2.7). Summing these gives 


(P.h = (Pt)j + (K)i 

(7.4.41) 

As in Chapter 6, a scheme which obeys the constraints 


(P+)j > min (-(Pr)i,0) 

(7.4.42) 

(Pf )j > min (-(P, + )j, o) 

(7.4.43) 


is second-order accurate. The only remaining difficulty is the selection of Qj+i which 
satisfies these constraints. 

For a particular cell, the value of (P+)y required to satisfy (7.4.42) may be greater 
than (P 4 + )j. The difference is defined as 

A(P+ )j = max (-(P," )j,(A + )i] - (P, + )j (7.4.44) 

Using (7.4.38) and an assumption of linearity, the change in can be computed. 

< 7 - 4 - 45 ) 

It is well known that smallest changes in 4>j+i (i n the sense of length) occur when A<f>j + i 
is chosen along the gradient. That is 

A ^+ i = ( 7 - 4 - 4 «) 

2 d< Pj+± 
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where a is a scalar quantity. When this choice is substituted into (7.4.45) the result is 

9 W)i 

(7.4.47) 


• _L 9< ^i + 1 

9(P£)i 




This value is used to compute a new value of 4>j+\- From <f>j+x, a new value of t> J+ i is 
computed using equation (7.4.13). A new value of 1 is computed using (7.4.22). The 
value of Vj + i is then used to compute a new value of fy+i using (7.2.11) through (7.2.15). 
Recall that the matrix Y' 1 is given in Chapter 2 and is evaluated at 9j+i* 


There will also be a change in (Pr)j+ l when g J+ i is updated. This change can be 

d( P+ ) • 

estimated using (7.4.39). It has already been shown that the two gradients j ; - — - and 


^i+i 

have the same sign, term by term. Consequently we can conclude that their 


9(Prh + 1 

Wj+i 

dot product is positive, their included angle less than 90 degrees. Thus, since the change 

in 4>i+i is along ^ , it has a positive component along ^ J — . This implies an 

1 9 < Pj+i 9 < Pj+\ 

increase in (P“)j+i. By implication, because the same procedure is carried out at each 

cell, (P“)j can be expected not to decrease. 


Due to the convexity of the entropy function, the actual change in (P+)j will be 
less than the linear approximation would suggest. Exceptions occur when roundoff error 
becomes significant. Since the change falls short of what is needed, iteration is required. 
The function (P+ )j is reevaluated, new gradients are computed and the definition of v 
changes. A subtlety is that <f>j+i changes too, a consequence of the changing definition of 
v. Iteration is continued until A is less than some tolerance. 

7.5 The Effect of Area Ratio and Metric Terms 


In the quasi-one-dimensional Euler equations (which describe flow through a nozzle), 
the effect of varying cross sectional area appears through a source term in the momen- 
tum equation. This corresponds to a pressure acting normal to the nozzle surface. In 
addition, the area ratio scales the usual fluxes in a manner analogous to metric terms 
in two-dimensional problems. In this section, the effect of the source term on the cross 
sectional area is investigated. 

A typical differencing scheme for solving the quasi-one-dimensional Euler equations is 


(q,t)j + y [{(<*/);+ A - («/)j- A } - hj [a i+ ^ 



(7.5.1) 
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where (a/) J+ ± = a J+ ± f- + 1 . The quantity a J+ 1 is the cross-sectional area of the nozzle at 
x j + 1 . The quantity Vj is the volume of the j th cell, typically Vj = aj( x j+i — x j-\ )• The 
vector h is given by 

(7.5.2) 

which is the source term for these equations. Notice that the details of how fj+i is com- 
puted are not given; the scheme is unspecified except for conservation and the treatment 
of the metric terms. 



The entropy inequality which corresponds to (7.5.1) is 


> 0 


Using the chain rule on the first term gives 

(«„««,.)/ + XT [(«F) J+ j - i ] > 0 


(7.5.3) 


(7.5.4) 


The factor (g,t)j can be replaced using equation (7.5.1) to give (after multiplying through 
by the volume) 

-( s >q)j [{W)j+± - ( a /)j-i} - h 3 {°j+i - a j-±}] + [( aF )}+± - ( aF )i-i] ^ 0 

(7.5.5) 

At this point it helps to note the following useful identity which may be verified by con- 
struction (perhaps using MACSYMA). 

S, q h = S, q f - F (7.5.6) 

This may be applied directly to equation (7.5.5) to yield 

“ F>] ( a 3+\ ~ °J-|) 

((*/),+ J - («/),. i) + UaF) j+ j - (aF)^,) > 0 

(7.5.7) 

Applying the definitions of af and aF and grouping terms gives the final form. 

«i + J - fi) + [F i+ i - Fj)] 

+ “j-i [-(S^it/i-Zj-jJ + lFy-Fj.,)] >0 (7.5.8) 


This form is independent of how the fluxes are computed and how the cell volume is 
computed. The source term does not appear. This agrees with classical thermodynamics; 
the source term is a pure force and therefore does not contribute to the entropy production 
rate. 
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7.6 Results and Discussion 


For the test problem, a converging, diverging nozzle geometry was used. The geometry 
is described by the formula 


a(x) = 0.8 * 



(7.6.1) 



Figure 7.1. — Nozzle geometry used for a test case. The con- 
ditions used had subsonic inflow and outflow with a supersonic 
region between x = 0.5 and x = 0.8. 


The particular case chosen has a steady solution with subsonic inflow and outflow. 
The flow is supersonic between the throat ( x = 0.5) and a shock which is located at x = 0.8. 

The first test was to verify that the semi-discrete scheme described in Section 7.4 
worked correctly. The steady state exact solution was used to test the scheme. The 
effectiveness of the gradient following technique can be seen in figure 7.2. The solid line 
represents ( P t )j , which uses linear averages of the entropy variables to compute 
This results in a negative entropy production rate over much of the domain. In order to 
use the stability conjecture of Chapter 2, a positive entropy production rate for each cell 
was sought. The diamonds represent the minimum entropy production rate guaranteed 
by constraints (7.4.42) and (7.4.43). The dashed line, which represents ( P,)j , shows what 
was actually achieved when (7.4.42) and (7.4.43) were met. 

The quantity q,t, computed by the semi-discrete scheme, found its largest values at 
the sonic point (x = 0.6) and at the shock (x = 0.8). The entropy production rate at the 
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shock was four orders of magnitude higher than any other point in the domain. This is to 
be expected from the physics. The near-zero entropy production rate at the sonic point is 
an artifact of the method; at sonic points, the constraints (7.4.20) and (7.4.21) translate 
to (P+), = (Prh = (P.)i = 0. 



Figure 7.2. — Semi-discrete entropy production rates for the 
exact solution. The solid line comes from using Vj + x = v. 

The chaindash line comes from modifying Vj + 1 to satisfy a cell 
entropy inequality. 


In Chapter 4, the question of time advance was explored. It was pointed out for 
example that explicit Euler is not a suitable time advance scheme for convective problems 
such as this one. The fully discrete entropy production rate for this scheme will fall short 
of the semi-discrete entropy production rate by an amount that is proportional to the 
change in the solution. This in turn is proportional to the time step. This conclusion was 
tested numerically, again using the exact solution for test data. In figure 7.3, the difference 
between the semi-discrete and fully discrete production rates is shown for various values 
of the time step. Notice that the loss in entropy production rate is considerably greater 
than the rate itself (fig. 7.2), especially around the sonic point. This would imply negative 
entropy production rates for the fully discrete scheme. The almost perfect proportionality 
of the curves (note the log scale) indicates that the uncertainty in the evaluation time 
of S, qq (4.2.8) makes little difference. It might as well have been evaluated at q n . The 
flattening of the curves at very low values is due to roundoff errors. 

In this work, implicit Euler was employed. In contrast to explicit Euler, it has the 
property that the fully discrete entropy production rate will always be greater than the 
semi-discrete entropy production rate. In principle this implies unconditional, nonlinear 
stability. For example, in figure 7.4, the entropy production due to time advance is shown 
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Figure 7.3. — Explicit Euler time advance destroys entropy. 

As Chapter 4 pointed out, the amount is roughly proportional 
to At. 

as a function of the time step. The entropy production rate increases as the time step 
increases but not proportionately. This is due to differences in the evaluation state of S, qq 
(4.3.8). This state is approximately q n+1 , which is different for each of the curves. 



Figure 7.4. — Implicit Euler time advance creates entropy. 
This effect is especially pronounced at larger time steps. 


In practice, there are several implementation difficulties, not all of which were dis- 
cussed in Section 4.6. One difficulty is that to obtain all the benefits of an implicit scheme 
requires solution of a system of algebraic equations that become increasingly nonlinear at 
larger values of At. Essentially this requires evaluation of the fluxes at a future time. Since 
any value of At can be used, this is equivalent to knowing the time-dependent solution 




in advance. Such knowledge is not generally available; indeed it would obviate the need 
for numerical methods. As far as nonlinear stability is concerned, implicit methods are 
restricted to time steps sufficiently small that linear approximations in q are valid. For 
example, the smallest time steps used in figure 7.4 were constrained by roundoff error, 
while the largest time steps were limited by the ability to converge the iterative process 
by which q n+1 was obtained. 


In practice, terms involving q n+1 are linearized about the current solution using a 
Taylor series. For example 


/w+1 




Vj ± i 

dqi 


(«? 





(7.4.7) 


The sum is taken over each unknown in the whole mesh, but many of the terms vanish. 
For this scheme, depends only on qj~i , qj,qj+i , <lj+ 2 - Since analytical expressions for 
the derivatives were difficult to obtain (to say the least), it was decided to evaluate the 
derivatives in (7.4.7) numerically by perturbing q slightly and measuring the the changes 
in / which result. This requires twelve additional evaluations of 1 per step (the number 
of variables per node times the width of the stencil). 


Once the derivatives are computed, a standard penta-diagonal matrix solver is used 
to advance in time. This gives an approximate value for q n+1 . For computing figure 7.4, 
this process was iterated until a very accurate value of g n+1 was obtained. On the other 
hand, for the main computation, the modification discussed in Section 4.5 was used to 
allow an accurate measure of the fully discrete entropy production rate. The modification 
has no effect at all when / is a linear function of q. 


It was noted with dismay that linearization in time has a disastrous effect on stability. 
For small time steps, the linearization is a good approximation to the nonlinear reality 
and the predicted behavior occurs. For Courant numbers greater than one, however, the 
linearization errors may be sufficient to cause local violations of the second law. These 
errors were largest during transients, especially those with moving shocks. As convergence 
approached, it became possible to take slightly larger time steps. With exact updates, in 
principle, there is no linearization error and no time step limitation. 

The last test consisted of starting with a fairly arbitrary initial condition and allowing 
it to evolve in time to reach a steady state. The initial condition chosen was obtained by 
linear interpolation of the conserved variables between the endpoints. A constant time step 
of 0.05 was used, roughly equivalent to a CFL number of |. The initial entropy production 
rate distribution is shown in figure 7.5. 

The solid line shows the fully discrete entropy production rate. Notice that it is 
greater than the semi-discrete entropy production rate (represented by the diamonds) 
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Figure 7.5. — Initial entropy production rate. The initial 
conditions are a simple linear interpolation of the boundary 
conditions. 

at each grid point. The difference is proportional to the change in q in the first time 
step. There is a slight glitch in the semi-discrete entropy production rate at x = 0.8, 
the eventual location of the shock. This glitch does not reflect numerical prescience, but 
rather a very slight metric discontinuity, an artifact of the mesh generation scheme (the 
mesh was almost equally spaced, but one point was placed at the shock location). This 
tiny glitch doesn’t appear in the fully discrete entropy production rate; it is swamped by 
the temporal dissipation (note the three order of magnitude difference in the two rates). 
Similar glitches at the boundaries are simply initial transients. To avoid roundoff problems, 
a small amount of entropy is intentionally produced in each cell. This leads to an artificial 
floor of 10' 10 on the entropy production rate. 

The exact solution was used to specify the conservative variables at inflow and outflow. 
This formally overspecifies the boundary conditions, so minor glitches at the boundaries 
can be expected. 

The initial condition was advanced to a nondimensional time of 25, at which time the 
solution appeared to be converged. The plot of density is shown in figure 7.6. The exact 
solution is shown as a solid line, and the computed solution is shown as the diamonds. 
Also plotted are the values of 1 from which the flux is computed (vertical ticks). Notice 
that none of these appear “in” the shock. 

Although the energy appears much the same as the density (and for that reason isn’t 
shown here) the momentum looks substantially different. It is shown in figure 7.7. Notice 
the substantial oscillation in the computed solution around the shock location. Although 
the computed solution oscillates, the values of </ J+ i are very close to the exact solution. 
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Figure 7.6. — Final density distribution. 



Figure 7.7. — Final momentum distribution. The computed 
solution matches the exact solution well on cell boundaries, 
where fluxes are balanced. It matches less well at cell centers, 
especially at the shock. Even though an interpolation is mono- 
tone in the entropy variables, it may not be monotone in the 
conservative variables. 

This brings up an interesting question concerning uniqueness of the solution. The 
solution is considered converged when it stops changing, i.e. when the fluxes reach their 
exact values. It seems likely that more than one solution is capable of interpolating the 
same values of 1 . Thus, there may be more than one solution. 
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The interpolation is limited by the consistency condition, which is applied to a partic- 
ular set of characteristic variables. Enforcing an interval boundedness condition in these 
variables does not imply such a condition in other variables, due to the nonlinearity of the 
transformation. The solution shown in figure 7.7 is a case in point; there is an interval 
immediately to the left of x = 0.8 for which ( pua)j + i does not fall between ( pua)j and 
(pua) j+1 . 

Other researchers 23 have suggested that one cannot expect interval boundedness in 
momentum at a shock on a finite mesh. The uncertainty in the shock location (it is only 
known to be somewhere in the interval) permits the computed state within the shock to 
vary within some limits. Its steady state value is a function of the initial conditions. The 
interval bounded quantities are p , it, and p. Under these conditions, interval boundedness 
for momentum occurs only in the degenerate case where the shock is reduced to a single 
interval. Ultimately, the problem stems from the approximation that qj is equal to the 
cell average. All schemes which use this approximation (including Godunov’s scheme and 
others) generate a similar glitch in the solution. 



Figure 7.8. — Final entropy production rate. Notice the 

strong peak at the shock (x = 0.8), as required by physics. The 
dip at the sonic point (x = 0.5) is an artifact of the numerical 
scheme. The low amplitude oscillations at the left boundary 
may be due to the overspecified boundary conditions. The 
agreement between semi-discrete and fully discrete values in- 
dicates very small changes in q. 


The final entropy production rate is shown in figure 7.8. Notice that the semi-discrete 
and fully discrete production rates correspond, as they must at convergence. Notice 
the four-order-of-magnitude spike in entropy production rate which occurs at the shock 

2 T. Barth, personal communications 

3 T. Barth, Some Notes on Shock Resolving Flux Functions, 1988 (unpublished) 


92 




(x = 0.8). The analytic solution shows a Dirac spike at this location, accounting for all the 
entropy produced in the domain. There is no oscillation in the computed entropy produc- 
tion rate near the shock to correspond with the observed oscillation in momentum. This 
is because steady state entropy production rates depend only on <fy+i values (the tiks in 
figure 7.7) which exhibit no such oscillations. Another feature of interest is a dramatic dip 
in entropy production around the sonic point (x = 0.5), where the constraints (7.4.20) and 
(7.4.21) evaluate to (P+ )j = ( P~)j = 0. Finally, some ringing in the entropy production 
rate is apparent at the inflow boundary. This is probably due to a slight inconsistency 
in the boundary conditions, which are overspecified. There is no apparent effect on the 
solution from this ringing. 

Overall then, it appears that satisfaction of a cell entropy inequality is sufficient to 
produce a stable scheme. Such a scheme has very little dissipation. On the bad side, it 
evidently isn’t sufficient to avoid the unphysical oscillations in momentum and it has a 
severe time step limitation. 

Despite its practical limitations, this scheme is interesting for its theoretical properties. 
Its computation of actual values for 1 is shared only by Godunov’s scheme, a first-order 
explicit scheme. The second law allows difficulties with grid spacing or time step to be 
pinpointed. This approach also provides insight to grid embedding, multigrid, and other 
interpolation-related problems. 
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Chapter 8. THE TVD CONNECTION 


8.1 Introduction 

This chapter reviews the basic concept of Total Variation Diminishing (TVD) schemes. 
Where applicable, it explores the connections between the TVD approach and the satis- 
faction of a cell entropy inequality. Many of the TVD ideas shown here are due to Harten 
(ref. 11). They are also found in an excellent survey paper by Sweby (ref. 27). 

Harten observed that existing schemes for solving convection problems were unsat- 
isfactory with regard to dissipation. One choice was monotone schemes (ref. 28). These 
schemes are guaranteed not to produce any new extrema. This property is useful in pre- 
venting negative densities for example, in computations involving very strong shocks. Un- 
fortunately, all linear monotone schemes are at best first-order accurate spatially (ref. 29). 

The other major choice consisted of adding constant coefficient dissipative terms to 
the PDEs as suggested by Lax (ref. 1). Since these are not part of the original equations, 
the coefficient is made as small as possible. Straightforward, symmetric differencing of the 
dissipative terms leads to stable schemes with second-order accuracy. An unfortunate side 
effect is the introduction of new extrema, oscillations, especially in the neighborhood of 
shocks. 

Thus algorithm designers were faced with a choice between first order methods which 
are robust, and higher order methods which aren’t. In either case, the dissipation is linear, 
easy to apply, and easy to analyze using Fourier analysis. 

A number of ad hoc methods emerged to make the dissipation coefficients dependent 
on the local solution. For example, Jameson suggested that the dissipation should consist 
of a linear second difference term multiplied by a term which depends on the second differ- 
ence of pressure (ref. 30). This “sensed” the presence of shocks and locally increased the 
dissipation there to remove oscillations. Pulliam (ref. 31) showed a qualitative similarity 
between this approach and upwind differencing. 

Another approach was taken by Boris and Book in their Flux Corrected Transport 
(FCT) approach (ref. 10). A second-order accurate scheme was used to compute fluxes 
between cells. Where these would result in new extrema, a first-order, monotone flux was 
also computed and an interpolation between the two was used. The interpolation was 
chosen in such a way as to eliminate the undesirable new extrema. 

This strategy was formalized by Harten in his TVD schemes (ref. 11). Harten aban- 
doned the requirement of monotonicity in order to achieve second order accuracy. Instead, 
he substituted the requirement that the Toteil Variation (TV) of the solution must not 
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increase with time. The concept of total variation is formally and unambiguously defined 
only for scalar equations in one space dimension. Given a solution vector u defined on J 
grid points, TV is defined as 

J 

= ( 8 . 1 . 1 ) 

;= i 

A TVD scheme is one in which this quantity can be relied on not to grow in time. 
Generally, periodic boundary conditions are assumed in the analysis. Very recently there 
has been some work on finite domain TVD schemes (ref. 32), but this is not discussed 
here. 


Harten explained how to construct schemes which are TVD. In particular, if the 
difference scheme chosen can be expressed in the form 

„»+> = - Cj.iW - «?_,) + U >+! («? +1 - «?) (8.1.2) 

where Cj_ i and are data-dependent scalar coefficients, and if 

> 0 
d h , > 0 

C j+ ,+D j+ t< 1 (8.1.3) 

then the scheme is TVD. The proof of this assertion is subtle but short. It is given in 
Appendix A in some detail. Conditions (8.1.3) are sufficient, but not necessary; there may 
be TVD schemes which do not meet these conditions. 

The concept of TVD schemes has been extended to implicit methods, through the use 
of additional conditions. It extends to multidimensions on regular grids through the use 
of operator splitting. Finally it extends to certain systems of equations through the use 
of eigensystem decomposition. A reasonable summary can be found in the work of Yee 
(ref. 33). 

The goal of the TVD approach is to produce solutions which are aesthetically pleasing. 
This is not necessarily the same as solutions that are correct (see fig. 5.4). The second 
law of thermodynamics provides a way to tell the two apart. 

8.2 Scalar Wave Equation 

This section explores the application of TVD ideas to the scalar wave equation. It 
also spells out the relationship between the TVD conditions and the requirements for a 
cell entropy inequality. 
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Consider the scalar wave equation: 


u,t + cu, x = 0, c > 0 


( 8 . 2 . 1 ) 


This is typically approximated by a central difference with a dissipation term 


u 


= „» _ - uJ_j) + ir[o j+ ,(«? +1 - «J) - a y . } («7 - «?_,)] 


( 8 . 2 . 2 ) 


The parameter a J+ 1 is a dissipation coefficient which may or may not vary from point to 
point. As written, equation (8.2,2) describes a conservative scheme which is consistent as 
long as a is less than unity and stable when a is positive. The parameter v is the CFL 
number: v = ^ L . 

At this point, (8.2.2) must be manipulated to look more like (8.1.2) so that the TVD 
conditions can be applied. The notation is simplified if u” +1 and are expressed in 

terms of differences. 


Au j+ i = u 7 + 1 ~ u 7 
Au j-\ = u 7 - u 7-i 

With these definitions (8.2.2) becomes 

“; +1 ±v(Au j+ i + Auj_i) + 1 Au h i - aj_iAv.j_i) 


By combining terms, this becomes 

,n+l - 1 


u j = «? - 2 ^! - ) Au j+\ - 2 U ^ + a i-0 Au j-$ 
Sweby (ref. 27) defines = 1 — a j+i- With this definition, (8.2.5) becomes 

u] +1 = uj - |^ i+ iAu i+ i - ii/( 2 - 


(8.2.3) 


(8.2.4) 


(8.2.5) 


( 8 . 2 . 6 ) 


By inspection, the only constant value of <j> which fits the TVD conditions (8.1.2) is 

<j> = 0. This results in the familiar first-order upwind scheme. There is, however an elegant 

Au i 

trick which improves accuracy. Sweby defines the quantity rj = A * * , and eliminates 


A«j + i from the equation. 


u 


n+1 

j 




(8.2.7) 
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which also has the the form of ( 8 . 1 . 2 ). The scheme is TVD if 



(8.2.8) 

This is equivalent, for positive v to 


~ Tj 3 » V 

(8.2.9) 

For v < 1 , tills becomes 

4>i+i 

2 < 3 > <j>; 1 < 2 

~ r i 3 * - 

(8.2.10) 

This can be separated further into the inequalities 


cs 

VI 

+ 

VI 

0 

(8.2.11) 

0 < <2 

Tj 

(8.2.12) 


These inequalities require that ^4.1 = 0 if rj < 0 . For rj > 0 , the inequalities are 
plotted in figure 8 . 1 . As long as <f>j+i{rj) lies within the shaded region, both inequalities 
will hold and the scheme will enjoy the TVD property. To achieve second-order accuracy 
in the limit as Aar — * 0 requires only that <£ J+ i(l) = 1 and that 4 >j + ^ be a continuous 
function of rj. (Several second-order accurate schemes are shown in figure 5 . 3 ) 



Figure 8 . 1 . — Schemes in the shaded region satisfy the TVD 
conditions ( 8 . 1 . 3 ). See figure 5.3 for some typical schemes. 
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There is a connection between TVD ideas and the entropy analysis discussed in pre- 
vious sections. For the wave equation the connection is known precisely. In Chapter 5, the 
semi-discrete entropy analysis was shown for this problem (eq. (5.1.8)). What is needed 
here is a fully discrete analysis using explicit Euler time advance. This is derived in Section 
5.4 (5.4.6) and is restated here. 


A t(P.)f B = v 





(8.2.13) 


If v is positive (the case where v is negative is a trivial extension), then a sufficient 
condition for positive entropy production rates is (5.4.9), plotted in figure 5.2. This is 
restated here for convenience. 


u 


-1 < 


£±i 


— u 


- < 1 


„ u . , 

3 j~\ 

It remains to define and Uj_ i. An elegant choice is (6.2.5) 


1 + u 


=u j+ 2^+* Au >+i 


(8.2.14) 


(8.2.15) 


and an analogous definition for With these definitions and the definition of rj 

previously given, inequalities (8.2.14) simplify to 


2 — 1 + v 


(8.2.16) 


These inequalities (for u = |) are plotted as the shaded region of figure 8.2. From 
previous analysis, the region between the two parallel lines satisfies inequalities which 
meet the TVD conditions (8.2.10). The area to the right of <f> = 2 is generally excluded 
on consistency grounds. Thus the cell entropy inequality is more restrictive than the TVD 
conditions in this case. Schemes which are consistent and satisfy the second law are TVD 
schemes. The converse is not true. 

8.3 Burgers’ Equation 

In this section, a TVD scheme is derived for Burgers’ equation. A nonstandard 
approach is used to facilitate comparison with the scheme derived in Section 6.3. A more 
standard approach will be demonstrated in the next section. 

Consider a generic conservative scheme for solving Burgers’ equation. 

U T' = U] - -/?-{) (8.3.1) 
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Figure 8.2. — Schemes between the parallel lines satisfy the 
TVD constraints. Schemes in the shaded region satisfy a cell 
entropy inequality. Schemes where 4>j+i > 2 are generally not 

used. Thus the entropy inequality is more restrictive than the 
TVD constraints. 


For Burgers’ equation, this reduces to 


Factoring the difference of two squares gives 

(8.3.2) 

“” +1 = “? •?-»*&! — *-» > 
This looks like the wave equation except that the wave speed is 

(8.3.3) 

c i - j ( u ?+$ +u 7- 1) 

(8.3.4) 

Defining vj = gives 

u 7 +1 =u ?~ 

As before, the definitions of tt? , and u? x are 

3+2 1 t 

(8.3.5) 

u l+ \ = u i 

(8.3.6) 

u l-\ 

(8.3.7) 
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With these definitions (8.3.5) becomes 


u 


n+l _ 


U 1 - 2 U ^j+\ Au 3+l " 2 Uj ( 2 ~ 


(8.3.8) 


which is virtually the same as (8.2.5) except for the definition of Uj. Thus the analysis of 
the previous section applies. If uj is positive, the scheme is TVD if 


—2 < - 4>i_i < 2(1 - ^ 


V ) 


(8.3.9) 


Recall that in Section 6.4, (u J+ ± — Uj) was assumed to be much less than Uj (leading 
to (6.4.8)). Here, as there, that assumption results in the approximation that Uj ~ u j ^ • 
Under this approximation, the scheme defined in Section 6.3 is identical to the TVD 
scheme just derived. Thus the connection between schemes which diminish total variation 
and those which produce entropy depends on the approximation of locally constant wave 
speed. Section 6.5 discusses the merits and drawbacks of such an approximation. 

8.4 Roe’s Approach 

There is another elegant way of deriving TVD schemes for Burgers’ equation. This 
approach, due to Roe, has the advantage of avoiding the approximation for Uj. More 
importantly it is easily (though not rigorously) extended to systems of PDEs. 

The approach taken by Roe is to use the intermediate value theorem to compute ± . 
That is 

fj+i fj = (/>u)j-|-i { u j+i u j) (8.4.1) 

where (f, u )j+i is evaluated at some value of u between and Uj. The intermediate 

value theorem guarantees that such a point exists. In particular, for Burgers’ equation, it 
is the linear average of uj and uj+ j. 


Substituting the definition of /(«), the left-hand side of (8.4.1) factors. 

\ [( u i+i ) 2 - ( u ; 2 )] = i^j+i - uj) 


(8.4.2) 


Equation (8.4.1) holds iff ( f, u )j+i = \( u i + 1 + Uj). Since /, u = u, (/, u ) J+ i should be 
evaluated at the linear average of uj and uj+ j. 

All of this can be incorporated into a TVD scheme as follows. Consider a central 
difference scheme with dissipation added, analogous to (8.2.2). 


u 


n +i — ,, n _ 'j 

~ U i 2Ax [Ij+1 /j " lJ 

At 


+ Ax [^ ,w a i+j Au j+$ (f' u )j-\ a j-\ Au j -\ ] 


(8.4.3) 
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Manipulating the central difference term gives 


1 At 


tt " +I = «? - sis wr+> - m + «r - /”-.)] 

+ ^ [(/'« tj+J a j+j Au )+i - Ur. 

Applying (8.4.1) directly, gives 


(8.4.4) 


_ 1 1 _ 1 At 

«7 +1 — U 1 ~ n~r~ 

3 J 2 Ax 

At 


(8.4.5) 


At this point the definition )j+i is introduced. Collecting terms gives. 

«i +1 = *? - («;+* + + K+* - ^h+i Att j+i ( 8 - 4 - 6 ) 

Finally, the substitution is made to express a in terms of <f> as for the wave equation 


u ] +1 = U 1 ~ \( 2 ~ “ ^ i +^ i +1 Au i + 

At this point an appropriate definition of rj is 


(8.4.7) 


r _ -i-t A "j-i 
’ *i+j A «i+i 

This can be used to eliminate Vj+i from (8.4.7), leading to 


(8.4.8) 


u? +1 = u? — i/,-_ l 
3 3 3 j 


(8.4.9) 


Provided that and i/ J+ i are both positive, this is in the same form as (8.2.7), 

from which a TVD scheme was constructed. A similar expression holds if both Vj-i and 
i ^ J+ 1 are negative. At shocks and sonic points (rj < 0) it is common to drop to first-order 
accuracy (i.e. <f> = 0). Unlike the scheme of Section 8.3, it is possible to evaluate i/ j+ i 
exactly, without iteration. 

The difference between Burgers’ equation and the linear wave equation is the treat- 
ment of the nonconstant wave speed. Not surprisingly, this is also one of the main dif- 
ferences between the various TVD and entropy producing schemes discussed so far. The 
other is the way that sign changes in the wave speed are handled. 
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8.5 One- Dimensional Gasdynamics, Roe's Identity 


In the previous section, the intermediate value theorem was used to factor a flux 
difference into an intermediate wave speed times a difference in the dependent variables. 
This approach generalizes to systems of equations as well. That is, there may exist a value 
of q ^ + 1 such that 

/i+i ~ fj — {fiq)j+x(Qj + 1 — Qj) (8.5.1) 

where (f,q)j+x means f, q evaluated at q J+ i. As shown by Harten, Lax, and van Leer 
(ref. 34), such a value exists if an entropy pair (5, F ) exists. Still, it may not be easy 
to find. Fortunately Roe has provided us with for the Euler equations (ref. 8). Its 
components are computed as follows: 


where, 


_ Pj+ly/Pj + Pjy/Pj+l 
Pj+ 3 \/Pi + y/Pi+1 

(8.5.2) 

__ u j+ly/Pj + U jy/Pj+1 

j+i y/P] + VP*' 

(8.5.3) 

k t _ hj+iy/pj + hj^/pj +i 

■ 7+3 y/Pj + y/Pj+\ 

(8.5.4) 

+ 

£ c 
i 

i 

in 

(8.5.5) 


/ ' 2 ' / i j ^ 1 2 ' j i j « v / 

The discovery of this identity allows the construction of a TVD scheme for systems in 
an elegant way. First, the scheme is written as a central difference scheme with dissipation 
to be added later. 

Qj +1 =9? ~ ~ f?- 1) + dissipation (8.5.6) 

Next the flux difference is replaced with two differences 


«" + ‘ =«?-§; [(/"+. -/”) + (/”- fU )] + di®»P*‘Son 


(8.5.7) 

Roe’s identity is used to replace the flux differences. 

9j +1 = 9j ~ ^ [(/>? )j+ 1 (?”+i ~qj) + (/„ )j- 1 (?? “ 1 )] + dissipation (8.5.8) 

The quantities (/, 9 )j+i and (/„ )j-i are matrices in this case. They can be replaced by 
the usual eigensystem decomposition 


«” + ‘ = «* - - »*> 

+ ,) + dissipation 


(8.5.9) 
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At this point, if all the eigensystems were the same, (i.e., = Yj) then 

the equation would formally diagonalize into three nonlinear, scalar wave equations which 
could be treated just like Burgers’ equation. In fact, this is not the case. This leads to an 
arbitrary choice in the pseudodiagonalization of the equations. Roe’s choice is to treat the 
quantity Y ml j + i(qJ +1 — q") as a difference in a new variable q. The equivalent equation is 

wr 1 = wj - -«) + a m(* - «-}>) ( 8 - 5 - io > 


The TVD techniques described in Sections 8.3 and 8.4 can then be used to find q J+ ^ 
and These are then multiplied by Yj + ± and Yj_^ respectively to recover f* + x an ^ 

fV_ i . Finally, q is advanced in time using 

3 i 


9 ; +1 - 




(8.5.11) 


Other, similar schemes are possible. All require ignoring the difference between eigen- 
vectors at nearby points. While numerical schemes for Burgers’ equation allowed some 
ambiguity in the speed of convection, schemes for the Euler equations aren’t even clear on 
precisely what is being convected. 

The fundamental problem goes deeper than approximation inaccuracies. The TVD 
property doesn’t hold for systems because of the interactions between the different eigen- 
vectors. This is true for the PDE as well as for the numerical schemes discussed here. 

8.6 Summary of the TVD Approach 

There are some similarities between the TVD approach and the entropy approach of 
the first seven chapters. For the linear wave equation, Section 8.2 showed that certain TVD 
schemes also satisfy a cell entropy inequality. The PDE in fact possesses both properties. 
On the other hand, in Section 5.5 it was shown that at least one TVD scheme does not 
satisfy even a global entropy inequality, leading to unsatisfactory results. Overall it would 
appear that the second law of thermodynamics is the fundamental property and that the 
TVD property is simply an important side effect. 

For Burgers’ equation, Sections 8.3 and 8.4 showed that it is possible to construct 
schemes which rigorously satisfy the TVD conditions. In Section 6.3, a scheme was con- 
structed which rigorously satisfied the cell entropy inequality, though this requires iteration 
at each cell and is very expensive. It is not clear whether such a scheme is rigorously TVD. 
It is clear from figure 6.7 that the particular TVD scheme used in Section 6.5 does not 
always satisfy a cell entropy inequality. This results in a shock that is slightly too strong 
at certain times, but the solution is otherwise quite satisfactory. Without the exact so- 
lution to compare with, the error would not have been detected. The main ambiguity in 
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nonlinear scalar equations of this type is the computation of a local wave speed and the 
variation of this wave speed over the cell. Indeed, it was shown that the assumption of a 
cellwise constant wave speed leads to the entropy errors just described. 

For one-dimensional gasdynamics, it is not possible to construct consistent schemes 
which formally diminish total variation. This is because the physical problem does not 
have the TVD property. Many schemes have been constructed using TVD ideas by ig- 
noring the interaction between the eigenvectors. These schemes appear to give excellent 
results, especially for reasonable ( M < 5) Mach numbers, in spite of the approximations in 
their derivation. For hypersonic mach numbers (Af > 5) the approximation of locally con- 
stant eigenstructure can become so bad that serious accuracy and convergence problems 
develop 1 . 

Although this section only deals with explicit TVD schemes, there has been consid- 
erable work done on implicit TVD schemes by Yee, Warming, and Harten (ref. 35). Such 
schemes produce results which are independent of the time step. Although the lineariza- 
tion is only an approximation to the gradient, it does allow larger time steps than the 
equivalent explicit method. 

A further modification avoids the need to explicitly calculate the eigenvectors at each 
grid point. This is a considerable savings in computational work at the expense of slightly 
more dissipation. Known as symmetric TVD, it was originally proposed by Davis (ref. 36), 
and developed further by Yee(ref. 37), and others. 

In two or three dimensions, the concept of total variation is not uniquely defined. 
Under one definition, the Lj norm of the solution should diminish. It has been shown 
(ref. 38) that schemes which diminish this norm in more than one dimension cannot be 
more than first-order accurate in space. The usual fix is to split dimensionally, applying 
a one-dimensional TVD operator along each of the coordinate directions. Such schemes 
are second-order accurate, but are not formally TVD in multidimensions, even for scalars. 
Furthermore, they are not easily applied to unstructured, finite element style grids. 

In problems with viscosity, heat transfer, or chemistry, there is no formal eigenvector 
decomposition as there is for the Euler equations. This makes it difficult to diagonalize the 
equations (even approximately) and apply the TVD property. The usual attack is to treat 
the inviscid convection terms separately and apply the TVD operator to them. The other 
terms are differenced in a symmetric stencil of suitable accuracy. The argument goes that 
such terms are inherently dissipative and hence can be relied upon not to create wiggles. 

The problem of extending the cell entropy inequality ideas to multi dimensions, un- 
structured grids, and more complete equation sets will be addressed in Chapter 10. There 
it will be shown that the second law extends nicely to all these cases. 

1 T. Barth, Some Notes on Shock Resolving Flux Functions, 1988 (unpublished) 
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Chapter 9. VARIOUS OTHER SCHEMES 


9.1 Flux-Based Schemes 

With the exception of the previous chapter, all schemes discussed so far have been 
variable-based schemes; i.e., schemes in which the numerical flux is determined by first 
computing the state of the independent variables at the midpoint. That is 

/j+3 = (9.1.1) 

This approach allows the computation of an entropy flux, that is consistent with the 
conservative flux, With the exception of Godunov’s scheme and its variants, there 

are no well-known schemes which use a flux of this type. 

It is more common to construct a numerical flux as a combination of the fluxes com- 
puted at the grid points. That is 


fi+i = M/,+1, fi) (9.1.2) 

A typical example is fj+\ ~ Wi + /i+i). the central difference flux. Such schemes are 
referred to here as flux-based schemes. 

Flux-based schemes are difficult to analyze from the point of view of a cell entropy 
inequality. This is because the entropy flux F is not a function of /. At a steady shock, 
for example, the conservative fluxes are constant through the shock, while the entropy flux 
jumps discontinuously. 

Because F is not a function of /, it is difficult to compute , which is needed to 
compute the entropy production rate for the cell. If entropy is considered at all, the usual 
approach is to independently approximate F. Although the global entropy production 
rate is unaffected by such an approximation, individual cells may have considerable errors. 
Such an approach is discussed in Section 9.3. 

There is another way to analyze the entropy production of flux-based schemes, a way 
which involves changing the quadrature rules by which fj+i is computed. Recall from 
Chapter 3 that the underlying integral equations for a convection problem are 

I q, t dv+<f /*n<L4 = 0 (9.1.3) 

Jv Jav 

Up to this point, one-dimensional problems have been treated as though the state were 
piecewise constant in each cell. In order to treat q as a continuous function, a profile has 
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Figure 9.1. — The volume integral behaves as though q(x) is 
piecewise constant over each control volume. Yet, the surface 
integrals are well defined. 


been introduced between cells (see fig. 9.1). This profile is assumed to be sufficiently 
narrow that the volume integral is unaffected, that is 



q, t dv = (q, t )j ;Ax 


(9.1.4) 


The surface integral however, is affected. The details of the profile determine the value of 
the state on the cell boundary. The surface integral is then evaluated as 


<t f >ndA = f(q j+ i)~ f{q ,_i) 
Jev 3 3 


(9.1.5) 


I 
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Figure 9.2. — In this view of a quasi-one-dimensional prob- 
lem, q(x,y) is piecewise constant over a slightly different area 
than the control volume. The thin, solid, curved lines are con- 
tours of q. The surface integrals (carried out on the dotted 
lines) must be evaluated in two parts. 
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It is often helpful to view one-dimensional processes as a limit of two-dimensional 
processes. For example, consider the convection in a two-dimensional pipe of unit cross 
section (fig. 9.2). If q = q(x ), q, y = 0, then a one-dimensional flow is described. Now, 
instead of allowing a profile to connect the states between cells, suppose the requirement 
that q, y — 0 is relaxed. Let the state be constant along a slightly different coordinate line, 
q, i.e., q, v = 0. This is shown as the thin solid line in figure 9.2. Assume that this line is 
sufficiently close to the cell boundary that the volume integral remains unaffected and can 
be evaluated as before. 

The surface integral contains one term for each cell face, as before. 

f av ^ = fj+\ ~ fj-\ (9.1.6) 

In 9.1.5, fj+i was evaluated easily because qj+i was a constant along the cell face. The 
situation shown in figure 9.2 is different, requiring the integration along each cell face to 
be done in two parts. 

fj+\ = \ [( 2 “ ^>-4 )/(?;) + ^+i/(9j+i)] 

= \ [( 2 “ ( t > H\ )fi + (9.1.7) 

The parameter depends on the geometrical details shown in figure 9.2, details which 
are adjustable. The flux / J+ 1 varies linearly between fj and fj+i as 4>j+i goes from 0 to 
2. When <t>j+i = 1, central differencing is recovered. 

The advantage of the physical point of view shown in figure 9.2 is that the entropy flux 
can now be computed. It is integrated along the same path as the other fluxes, resulting 
in the expression 

F j+\ = \ [( 2 ~ 4>j+\) F i + (9.1.8) 

where F i+i = F(qj+i), etc. This allows the cell entropy production rate to be computed. 
It is sometimes possible to select in such a way that this (P t )j > 0. 

For an example, consider the wave equation again. What is the entropy production 
rate, given this new definition of FI (Since f — cu in this case, the new definition of / has 
no effect.) Following the same derivation path a s before (Section 5.1), the semi-discrete 
entropy production rate is 

Ax ( F »)j = ~ c “ u j) 2 - ( 2 “ ~ u i- 1) 2 ] > 0 (9.1.9) 

As before, the substitution rj = is helpful. For positive wave speeds this leads to 

(9.1.10) 
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Variable Based 

Flux Based 



Figure 9.3. — Both limiters are second-order accurate. Both 
satisfy an entropy inequality. The variable-based limiter has 
about half the dissipation of the flux-based limiter. 


which is satisfied if 

h xnin(r,- a ,l) (9.1.11) 

For comparison, the equivalent expression for the variable-based scheme is 

= min(|rj|,l) (9.1.12) 

These two limiters are plotted in figure 9.3. Both limiters are asymptotically second- 
order accurate because both are continuous through the point (1,1) on the graph. However, 
the flux-based limiter is about twice as dissipative as the variable-based limiter. 

For the wave equation with c = 1, fj+i = u j+±' Thus, every flux-based scheme can 
be viewed as a variable-based scheme and vice versa. Both limiters are stable and give 
satisfactory results. Thus, the difference between the two is purely one of analysis. 

Another problem arises when flux-based schemes are extended to systems of equations. 
The view shown in figure 9.2 allows only one interpolation parameter at each cell. This is 
too restrictive to guarantee positive entropy production rates at sonic points. Numerical 
experiments revealed significant problems at such points in the form of oscillations, glitches 
and instabilities. 

Flux-based schemes, as defined in this section, are an interesting dead-end. To make 
any further progress, it is necessary to somehow remove the restriction of having only one 
interpolation parameter per cell. While this is possible conceptually, rigorous analysis is 
difficult. 
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9.2 Flux Splitting and the Second Law 

An interesting connection exists between flux splitting and the second law of thermo- 
dynamics. It involves a Taylor series approximation and is therefore inexact. Nevertheless, 
some significant insights are possible. 

Formally, flux splitting only applies to certain purely hyperbolic systems of equations 
such as the Euler equations. These are written in the usual way 

+ /,. = 0 (9.2.1) 

The semi-discrete equations can be written 

4 ,. + -4-1) = o (9.2.2) 


The general idea behind flux splitting is to use the homogeneous property (fj = 
(f, q )jqj) or Roe’s identity (8.5.1), and the known eigensystem of f, q (7.3.2), to approximate 
a system of equations by several scalar equations. Methods which use fj = ( fiq)jQj are 
known as flux vector splitting methods while those which use Roe’s identity are known as 
flux difference splitting methods. 


As previously derived (3.4.3), the semi-discrete cell entropy inequality is equivalent to 


[-(s»M/> +} - fi) + (*■/+} - *») + HSMfi - f}-0 + (*> - *!-})] > o (9.2.3) 


The quantity on the left is the entropy production rate for the j th cell. In smooth regions 
it should approximately vanish. This notion can be formalized by expanding in a Taylor 
series about qj . For example, 


/j+i w fj + (/>«); (fy+i — 9j) + o (?i+$ 9j) T (/w«)i(9j+i 9j) 


(9.2.4) 


The term (qj+± — qj) appears so often in expansions of this type that it is replaced 
with the letter a. Similarly, b = (qj — fy_i). Finally, the subscript j, now redundant, is 
dropped. Grouping the terms by powers of a and 6, (9.2.3) becomes 

-[SM-M + iFj-F,)] 

+ [—Siqfiq( a ) + ^w( fl )3 

+ “ [— fiqqO,) + a Fiqqd] 

+ b T F m t) > 0 (9.2.5) 
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The first two terms vanish by inspection. The second two terms vanish because of the 
compatability condition, F, q = S, q f, q . Cancelling the factor of (9.2.5) becomes 


[a T (F„, - S, ,/„,)«] - [i T (F,„, - S„f,„)b] > 0 


(9.2.6) 


This simplifies substantially since 


Fiqq — i^^qf w)>9 — &,qqf,q + S,qf 


iqq 


(9.2.7) 


which comes from the compatability condition and the chain rule. This cancels the f, qq 
term (a third-rank tensor) and (9.2.6) collapses to 


a T S,qqf,q<l b^ S \qq f \qb > 0 


(9.2.8) 


Since S, qq is negative definite, it is clear that f, q plays a key role in the satisfaction of the 
second law. This observation is clarified if the identity (7.3.10) is applied. Recall, this is 

s,„f„ = ~(Y-') T AY-' (9.2.9) 

Inequality (9.2.6) can be formally diagonalized using this identity. This allows the two 
terms to be combined and leads to the expression 

-*i) 

+ A2(5.2 — ^2 ) 

+A 3 (a| - b 2 3 ) > 0 (9.2.10) 


Satisfaction of (9.2.10) is a fairly simple matter, accomplished on a term-by-term 
basis. Thus, a system of equations is broken into three, nonlinear, scalar equations. The 
resulting scheme is identical to flux vector splitting (ref. 9). This allows the analysis of 
flux vector splitting from the point of view of the second law. 

In Chapter 6, the behavior of nonlinear scalar equations was studied. In particular, 
figure 6.2 show that the approximations Uj+i = Uj and = Uj + 1 may both violate 

a cell entropy inequality at sonic points. First-order flux vector splitting is constrained 
to use one of these approximations. It seems likely that a local violation of the second 
law is responsible for the sonic point glitch which is characteristic of first-order flux vector 
splitting. 

As a second-order scheme, will approach the average of Uj and Uj+i as the mesh 
is refined. Since the Taylor series is taken about uj , the Taylor Series errors do not vanish. 
Thus, the validity of this analysis for second-order flux vector splitting is questionable. 
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Flux difference splitting, on the other hand is generally linearized about a state that 
is, in some sense, the average of uj and v,j+ \ (perhaps the Roe average of Section 5.8). 
In addition to avoiding the sonic point problems experienced by flux vector splitting, the 
linearization error vanishes rapidly in the limit of mesh refinement. 

It is worth noting here that attempts were made (unsuccessfully) to construct a 
scheme which satisfied a cell entropy inequality using (9.2.10). Because of truncation 
errors, (9.2.10) does not imply satisfaction of (9.2.3). These discrepancies proved large 
enough to cause instabilities. 

9.3 Tadmor’s Entropy Scheme 

So far, the emphasis of this work has been on the satisfaction of the second law in 
each cell. On the other hand the only formal theorem involving entropy requires satisfying 
the second law over the whole domain. This section examines a scheme which satisfies a 
global entropy inequality without satisfying a local entropy inequality. 

Tadmor showed (ref. 15) that many equations of interest can be discretized in a way 
that globally conserves entropy, in addition to mass momentum and energy. He does 
not recommend such schemes but shows how to construct them. He also gives a fairly 
general form of dissipation which results in positive (global) entropy production rates. 
Such schemes fall into a class he terms entropy-stable schemes. 

In this section, the limiting case of entropy conservation is explored. Since Burgers’ 
equation produces entropy at shocks, an entropy conserving scheme is certain to be at 
variance with the physics. It is interesting to observe the behavior of such a scheme. 

The equation of interest here is Burgers’ equation. 

tM + ^(« 2 ),* = 0 (9.3.1) 

By carefully following the analysis in (ref. 15), the following flux is derived 

\ 

fj+i = ^(«j 2 + + Uj+i 2 ) (9.3.2) 

It is proved by Tadmor that this flux conserves entropy in the sense that 

(S„)i + i[FV +i - FV.j] = o (9.3.3) 

The quantity F*j+i (given in (ref. 15)) is an approximation to the entropy flux. It is 
given in general by the expression 

= jIFf+i - (S„W/, + > - f i+i ) ] + i[Fj + ~ A)1 (9.3.4) 
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which for Burgers’ equation reduces to 

F * j+ \ = + Uj) (9.3.5) 

Notice that this becomes exact, as does / J+ i when tij+i = Uj. Thus, both F *j+i 
and fj+i meet the Lax condition for consistency. Nevertheless they are independently 
approximated. That is, there is no 1 for which, simultaneously, fj+i = f( U 3+\) and 
F* j+i3 =F{u j+ ±). 

Using implicit Euler time advance to advance the semi-discrete scheme of (9.3.2) 
ensures that the fully discrete entropy production rate will be greater than the semi- 
discrete entropy production rate (see Section 4.3). The results for this scheme are shown 
in figures 9.4 and 9.5. Figure 9.4 shows the excellent agreement between Tadmor’s scheme 
and the analytical solution prior to the formation of a shock. On the other hand, the 
post-shock performance is terrible. Figure 9.5 shows the computed versus exact solution 
at the time of peak shock strength. Notice the severe oscillations in the neighborhood of 
the shock. 
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Figure 9.4. — Tadmor’s solution at onset of shock. 

In spite of the oscillations and generally poor quality, the solution remains bounded. 
In particular, the sum of all the tiy’s is conserved as is the sum of all the itj 2 ’s. Since 
5 = —u 2 for this problem the second law is satisfied globally. The question of cell-by-cell 
entropy production remains. 

Although (9.3.4) is a good approximation to the entropy flux in smooth regions, it is a 
poor approximation near the shock. It is also a poor approximation amongst the wiggles. 
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Figure 9.5. — Tadmor’s solution at maximum shock strength. 
Both u and u 2 are conserved quantities when, in fact, u 2 should 
be decreasing. Hence the wiggles. 


To show this, one may compute an exact by computing fj+i from (9.3.2) and then 
solving fj + x = |u J+ i 2 for u J+ i . The sign is given by consistency, that is has the 

same sign as uj + uj +i . Once the exact u J+ 1 is known, Fj+i can be computed using the 
formula 



(9.3.7) 



Figure 9.0. — Entropy production rate for figure 9.5. Notice 
that no net entropy is produced over the domain. There are 
individual cells that destroy entropy. These are the source of 
the wiggles in u. 
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With this definition for it is possible to evaluate the true entropy production 

rate, cell-by-cell, as in Chapter 6. This is shown in figure 9.6 for the solution shown in figure 
9.5. The total entropy production rate for the domain is zero, but the entropy production 
rate in some cells is negative. The entropy production rate computed by Tadmor is positive 
and very small, invisible on this scale, all of it due to the time advance. 

This discrepancy between two ways of computing entropy production rates is only 
evident in the presence of shocks. It is significant that (9.3.2) works quite well until a 
shock appears. At the onset of the shock, two things happen at once; oscillations appear 
in the solution and cells near the shock begin to display a negative entropy production 
rate. 


It must be noted that neither Tadmor or Osher ever intended for the scheme to be 
used in this way. In particular, Osher recommends the addition of a suitable limiter, 
to enforce a TVD condition. Along the same lines, Tadmor suggests the addition of a 
dissipative term of a particular form. This form is sufficiently general to allow the limiter 
recommended by Osher. Either way, the added dissipation increases the global entropy 
production rate. This (according to Tadmor (ref. 15)) allows the scheme to converge to 
the correct answer as the mesh is refined. 

The scheme is nevertheless instructive as written. It shows the one-to-one corre- 
spondence between local violations of the second law and unphysical features of the flow 
(oscillations). This in turn shows the importance of satisfying a local entropy condition in 
addition to a global one. 


9.4 Entropy and Wiggles 

The nonlinear stability conjecture of Chapter 2 depends on satisfaction of the second 
law, but only in a global sense. Throughout this work, the much more restrictive cell 
entropy inequality has been stressed. Certainly, if entropy is produced in each cell, it is 
produced in the domain as a whole. Thus, the satisfaction of a cell entropy inequality 
implies nonlinear stability, at least for scalars. 

It is definitely easier to construct schemes which only have to satisfy a global entropy 
inequality. For example it should be fairly simple to solve for fy+i such that (P+ )j = 
~(P.~ )j+i, since both of these depend only on q J+ 1 . Yet, it only makes sense to enforce 
the first and second law on the same control volumes, the computational cells. 

Perhaps the cell entropy inequality provides some desirable property in addition to 
stability, in exchange for the increased computational complexity. There are indications 
that this might be so. Recall that in Chapter 5, the schemes which satisfied a cell entropy 
inequality turned out to be TVD as well. The previous section showed a scheme which 
did not satisfy a cell entropy inequality but did satisfy the second law in a global sense. 
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The result was a scheme which was stable as promised, but with severe oscillations in the 
solution. This suggests that the cell entropy inequality provides control of oscillations, in 
addition to stability. 

Consider now a further illustration of this hypothesis. The problem considered was 
previously discussed in Section 7.6. The solution, given in figure 9.7, was obtained with a 
flux-based central difference scheme using a constant coefficient, second difference, dissi- 
pation model. The analytical solution for the density is given by the solid line while the 
computed solution is plotted as symbols. Notice the overshoots in the neighborhood of 
the shock. Intuitively these seem to be unphysical, yet the method is fully conservative of 
mass, momentum, and energy. Unfortunately, it does not satisfy the second law. 



Figure 9.7. — Compare this density distribution for a prim- 
itive one-dimensional Euler solution with that shown in figure 
7.6. Notice the oscillations in the computed solution around 
the shock. 


For steady solutions the entropy inequality simplifies to 

P t =A F • n <L4 > 0 (9.4.1) 

Jav 

where the quantity F is the local entropy flux. This allows an after-the-fact evaluation of 
solutions to the steady one-dimensional Euler equations with regard to their satisfaction 
of the second law. For this case the entropy flux F is (pua)s where s is the specific entropy 
of the working fluid and pua is the mass flow rate. The second law in this case amounts 
to a requirement of monotonicity of s, since pua is a constant for steady state solutions. 
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For steady solutions then, equation (9.4.1) should be obeyed over any control volume. 
The integrals reduce to differences in the one-dimensional case. By evaluating P t for 
control volumes bounded by adjacent grid points, it is possible to see whether this method 
satisfies a cell entropy inequality. In figure 9.8, the entropy flux is given. 
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Figure 9.8. — Plot of entropy flux for a primitive Euler solu- 
tion. The computed entropy flux is not a monotone increasing 
function of x. Therefore this solution does not satisfy the sec- 
ond law. Notice that the oscillations in entropy flux extend 
much farther than the oscillations in density. 
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As in figure 9.7, the symbols represent the computed solution. It is immediately 
apparent that the solution fails to satisfy the second law ( F declines) in approximately 
half the cells. It is possible to group some of the control volumes together to form new 
control volumes, each of which obeys the second law. Some grid points fall entirely within 
these new control volumes. Such points are considered numerical artifacts and are not 
plotted. It is possible to choose a (nonunique) subset of the data which does satisfy the 
steady state second law for each cell (i.e., one for which F is monotone increasing). Figure 
9.9 shows the entropy flux for such a subset. 

In figure 9.10, the corresponding density values are plotted for these points. The . 
wiggles are gone and the solution doesn’t look too bad, though the shock is smeared. 

In summary, this numerical scheme uses grid points inefficiently. The solution has 
roughly the same information content as a more dissipative scheme with half as many 
points. 

Incidentally, the choice of which half of the solution to use is not arbitrary as might 
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Figure 0.9. — Plot of entropy flux for selected grid points. 
These points form a set for which entropy flux is monotone 
increasing. The second law is satisfied over these points. 
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Figure 9.10. — Plot of density for selected grid points. This 
is the same set of points as in figure 9.9. Notice the absence of 
oscillations around the shock. 


be suggested by an odd-even decoupling argument. The points which are excluded from 
figure 9.10 are shown in figure 9.11. Although the solution appears reasonable, and has 
a sharp shock, it represents a physically impossible flow. The compression requires less 
work than the theoretical minimum and the post-shock expansion produces more work 
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Figure 9.11. — The density values excluded from figure 9.10. 

These values form a smooth set. However the shock is too 
strong and the minimum density is too low. This could lead to 
negative densities in practice. 

than the theoretical maximum. All this excess availability is consumed in the shock, which 
is noticeably stronger than it should be. As a practical matter, the minimum density is 
significantly low. When capturing a very strong shock, this could lead to negative densities 
or pressures. Very few schemes can tolerate this unphysical situation. 

Increasing the dissipation coefficient sufficiently will give a solution without any wig- 
gles on the original grid. Such a solution is shown in figure 9.12. When the entropy flux is 
plotted in figure 9.13, it becomes apparent that the entropy inequality is not satisfied over 
a large part of the domain. 

The trick used in the previous example would require removal of the entire supersonic 
region (indeed the accuracy is quite poor in that region). Thus, smooth solutions do 
not necessarily represent the physics accurately. They may contain less information than 
solutions with oscillations. 

At least for this example, the hypothesis appears to hold. If a solution satisfies, in 
each cell, conservation of mass momentum and energy and the second law, it will not 
exhibit unphysical oscillations. Such a principle is very useful in constructing numerical 
schemes. 
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Figure 9.12. — Density distribution using a highly dissipative 
scheme. This solution is smooth, but highly inaccurate. 
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Figure 9.13. — Plot of entropy flux corresponding to figure 
9.12. Although the solution in figure 9.12 is smooth, it does 
not satisfy the second law. Smoothness is not enough, even 
with conservation. 
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Chapter 10. PROPOSED EXTENSIONS AND SOME OBSERVATIONS 


The discipline of thermodynamics is widely applicable. Essentially it applies over the 
whole of continuum mechanics. This includes phenomena involving heat transfer, viscosity, 
and chemistry, in two or three dimensions. 

The numerical idea of a cell entropy inequality derives from the analytic ideas of 
thermodynamics. It seems reasonable that it should be just as widely applicable. This 
section discusses the details of how various extensions might be implemented. 

10.1 Multidimensions 


If the numerical techniques described in previous sections are to be useful, they must 
be extended to more than one dimension. This can, of course, be done in the usual 
way — by restricting the computation to regular grids and applying the one-dimensional 
algorithm in each dimension separately. This technique has widely recognized difficulties. 
Fortunately there is a more satisfying approach available. 


Recall from Chapter 3 that the particular equations of interest, in integral form, are 
the conservation law 

f q[t + At]dv- [ q[t]dv+ f <f { ndAdr = 0 (10.1.1) 

Jv Jv Jt Jsv 

and the corresponding statement of the second law 

»t+At ft t ft+^t f 

/ / P t dv dr = / S'tf + Afjdu- / S[t]dv + / <p ¥-ndAdr>0 (10.1.2) 

Jt Jv Jv Jv Jt Jav 


As a first step, assume that the domain of interest is tessellated into a finite number 
of control volumes, each of which will therefore have a finite size. Assuming quadrilateral 
cells arranged in a regular array as in figure 10.1, these can be identified by the two indices 
j and A:. The surface surrounding this control volume is also identified by the indices jk. 


The integral equation which applies to this volume is 

pt-fAt 


[ q[t + At]dv— f q[t]dv+ f f>ndL4dr = 0 (10.1.3) 

Jv jh Jv jh Jt JaVj h 

Because each control volume has a finite extent, the volume integral must be a continuous 
function of time. Dividing (10.1.3) by At and proceeding to the limit At — * 0 gives the 
semi-discrete form. 
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(10.1.4) 
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Figure 10.1. — A typical control volume in two dimensions. 


Equation (10.1.4) holds at each time. At this point, the equations can be simplified by a 
definition. 


1 




= ?jfc 


(10.1.5) 


With this definition, (10.1.4) becomes 


Vjk 


dqjh 

dt 


+ 


/ 

JBV jk 


f • n dA = 0 


( 10 . 1 . 6 ) 


In one dimension the surface integrals could be carried out exactly. This is not true 
in two dimensions. It is possible to approximate these integrals however. The control 
surface, dVjk, which surrounds the control volume, can be broken into four segments in 
a straightforward way. As labeled in the figure, q is assumed to be constant along each 
of the four segments. Referring to figure 10.1, these four values are qj+i Qj-i ki 

and qjk- 1- Subject to this assumption 

f • ndA = a j+ 1 k f j+ i k -n + a j _i fe f y _i k • n 

' n + a jk-\tjk-i -n (10.1.7) 

where in each term, n is an outward facing unit normal at the appropriate face. The cell 
face areas, e.g., a,j+i appear separately in each term, allowing for a more general grid 
than that shown. Combining (10.1.6) and (10.1.7), and dividing by the volume (also called 


/ 

Jb v ih 


121 



the metric Jacobian) gives 




n + 


dj-1 fc 


f; 


LXl ‘ 


1 fc # 

Vjk 3 3 

, a i‘-h 

n ^ f*L i • 


n 


( 10 . 1 . 8 ) 


For the rectangular mesh shown in figure 10.1, certain simplifications are possible. 
Using the traditional unit vectors i and j in the x and y directions, the dot product f • n 
can be carried out. Since the cell faces are aligned with coordinate directions, the following 
definitions are sufficient. 

/ = fi and </ = f • j (10.1.9) 

Since the various areas are Ax or Ay and the cell volume is Ax Ay, equation (10.1.8) 
can be written, for rectangular meshes, as 


( 9jk),t + 


fj+\k fj-\k t 9jk+\ 9jk-l _ 


Ax 


+ 


Ay 


( 10 . 1 . 10 ) 


This is analogous to equation (3.1.6), the one-dimensional version. The approxima- 
tion in (10.1.7) remains, (10.1.10) is not exact. This is the major difference between one 
dimension and multidimensions. 

The piecewise constant approximation used in one dimension is retained. 

q — qj k for all q within control volume jk (10.1.11) 

where qj k is the average value of q over the control volume. As shown in Section 2.2, this 
assumption maximizes the entropy within the control volume, given the known value of 
qjk- Overall then, two dimensions requires one more assumption (approximation) than one 
dimension. In addition to assuming q is piecewise constant within a control volume (as in 
one dimension), it is assumed to piecewise constant along the control surface. 

A semi-discrete version of the second law can be constructed in exactly the same way 
that (10.1.8) was constructed. This is 


(P.U = + ^F, +i k 


n + 


a j-i *_ 

_i_ a_ f, 


Vjk 




n 


a.* l i i a ■ l_ i 

-LaLLF ifc+ i • n + fc _ i • n > 0 

vjk 3 +3 Vjk 3k 3 


(10.1.12) 
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This does not involve any further assumptions. Since q is assumed piecewise constant 
along the cell boundary, P is computed as easily as f. For rectangular meshes this also 
simplifies. Using the definitions 


F = F • i and G = F • j 

inequality (10.1.12) can be written, for rectangular meshes, as 


(P»)jk ~ (Sjk),t + 


F i+\h-Fj-lh | Gjfc+j ~ > 

A* A y ~ 


0 


(10.1.13) 


(10.1.14) 


which is analogous to (3.1.8). The time derivative can be eliminated {is before by using 
the chain rule. In particular S,t = S, q T q, t , where q, t appears in (10.1.10). Making the 
substitutions yields 

(P.)ik = + i ' it+ ' Ay !7 ^ i ) 


~ . G i<-+\ Z G > t -< > n 

Ax Ay ~ 


(10.1.15) 


This is algebraically equivalent to 


(A ) » = 


Ax 


Ay [ 


- ( s jk),q T (fj+ik -fjk) + (-Fj+Ifc - Fjk) 
~ (Sjk)iq ( fjk — fj-ik) + (Fjk — Fj- Ifc) 

- (9jk + 1 — 9jk) + (Gjk+i ~ Gjk) 

- (S jk )„ T (g ik - , JM ) + (G it - G jk _ { ) 


> 0 


(10.1.16) 


which is similar in form to (3.4.3) but with the metric terms left in. This is more easily 
manipulated in the form 


(P»)jk - (F,)j+ife + (Pt)j-ik + (P*)jk+$ + (P»)jk-± - 0 

where, for example 

(A) i+ i* = A [-($*), - M + (F j+kk - F i4 )] 


(10.1.17) 


(10.1.18) 


Following the development in one dimension, a cell entropy inequality can be sat- 
isfied in several different ways. If the simplicity of a first-order scheme is desired, then 
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9j+i *>9jfc+i> and qj k _± must each be adjusted in such a way that the individ- 

ual terms of (10.1.17) are nonnegative. The usual approach of dimensioned splitting is 
equivalent to a pair of constraints. 

(A )/+!*+(>.);- 0 and (?,Wj + > « (10.1.19) 

The first constraint corresponds to the x direction and the second corresponds to the y 
direction. These can be satisfied using the one-dimensional techniques already shown. 
Such an approach will be second-order accurate in each of the coordinate directions, but 
may have serious anisotropic properties which are generally undesirable. 

In fact, it is only the sum ( P t )jk which is physically constrained to be nonnegative. 
This constraint can be met by using a two-dimensional analogue of the one-dimensional 
technique. First, some second-order accurate states such as qj+i * are chosen, as was done 
in one dimension. This can be done by linear averaging, perhaps in the entropy variables. 
This defines the bar terms such as (P,) J+ i fe . Summing these gives 

( 5 )j» = (a)i +1 * + (a) i - 1 *+(a)i» + 1 +(A)j»-j >o <10.1.20) 


If ( P,)jk > 0, then (10.1.17) is met by the average states and no adjustment is 
necessary. If ( P$)jk < 0 there must be one or more negative terms. Each can be reduced 
in magnitude to ensure ( P,)jk > 0. Let the sum of the negative terms be (i*)Jfc. This 
sum must be reduced in magnitude (multiplied) by some positive constant a, 0 < a < 1. 
This constant is given by _ 


(P.hk 

(P.)Jk 


( 10 . 1 . 21 ) 


The total change can be divided amongst the negative terms, with those of larger magni- 
tude getting a larger share of the change. For example, if (P,)j+Xk < 0 then 1 fc needs 
to be changed in such a way that 


( 10 . 1 . 22 ) 


where 



(10.1.23) 


and is the sum of the squares of the negative terms. The desired changes to interme- 

diate states such as q J+ i fc can be accomplished by gradient methods already introduced. 


An example suffices to show that this can be a significant improvement over dimen- 
sional splitting. Suppose that (P t )jk > 0, but (P a ) J+ i fc + (P«)j_i* < 0* I n this case no 
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dissipation needs to be added because (10.1.16) is already satisfied. A dimensionally split 
algorithm would add dissipation anyway because the first inequality of (10.1.19) is not 
satisfied. In fact, dimensional splitting always requires at least as much dissipation as the 
more general method of (10.1.22). 

10.2 Unstructured Grids 

There has been considerable interest recently in unstructured grids. Such grids are 
substantially easier to generate about complicated geometries than grids which are topo- 
logically rectangular. It is difficult to define a suitable dissipation operator on an unstruc- 
tured grid with existing methods. For example, dimensioned splitting is not an option so 
the usual definitions of total variation or fourth difference operators do not apply. 



Figure 10.2. — A typical control volume for an unstructured 
grid. 


The simplifications used for rectangular grids do not apply to unstructured grids, one 
cell of which is shown in figure 10.2. The control volumes are triangles in this case, but 
they could be quadrilaterals, pentagons or any other polygons. There is no requirement 
that they be convex or have the same number of sides, as long as they tessellate the 
domain. The distinguishing feature of unstructured meshes is that there is no obvious 
numbering scheme, as there is for logically rectangular meshes. This being the case the 
control volumes are usually identified by a single index, in this case j. Each control volume 
has a certain number of neighbors (three in this case). The control surface between the 
i th and j th control volume is denoted by the subscript ij. In figure 10.2, the neighbors of 
cell j are cells 1,2, and 3 for notations! convenience. 
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The derivation in Section 10.1 can be followed up to equation (10.1.6). At this point 
it is useful to note that 

/ f, ■ ndA = 0 (10.2.1) 

Jav, 

because the divergence of a constant vector vanishes. Adjusting for the difference in 
notation and subtracting (10.2.1), (10.1.8) can be written 

(?,),. + ^(fji - *» • n + ^(f ;j - fj) • n + - f, ) • n = 0 (10.2.2) 

Vj V j v j 

The dot product is carried out for each term, leading to the definitions 

Afji = aji(fji ~ fj) • n 

Afj 2 = aji({j2 — fj) • n 

Afj3 = aj3(fj3 — fj) • n (10.2.3) 

With these definitions (10.2.2) simplifies to 

kiU + y Mu + M* + Mjt) = o (10.2.4) 

In a completely analogous manner, the semi-discrete entropy inequality can be derived. 

{P.)j = (SjU + ^ [AF,-, + A F j2 + A F j3 ] > 0 (10.2.5) 

V 3 

As before, qj is assumed to be constant over the control volume and qij is assumed 
to be constant over each edge. Using the chain rule again (5, t = 5„q,t), it is possible 
to eliminate (5j),t using (10.2.2). This gives, for unstructured meshes, an expression 
corresponding to (10.1.16). 

v j(P»)j = — [Afji + A/j 2 + Afj 3 ] + [AFjj + AFj2 + AFj 3 ] > 0 (10.2.6) 

Combining these as before gives 

Vj(P.)j - [-(S„)l A /yi + A^i]+[-(S„)iA/ j2 + AFj 2 ]+|-(S„),A/jj + AF,-j] (10.2.7) 
Finally 

(P.)j = + (P,)j 2 + ( P.) j3 > 0 (10.2.8) 

where, for example 

(P.)yi s i [-(SJjA/fl + AFji] (10.2.9) 

The remainder of the derivation is identical to that of Section 10.1, except that the number 
of terms in this case is three instead of four. The intermediate states qj\ , qj 2 , and qj 3 can 
be constructed so as to satisfy (10.2.8). 

It is possible to recover the expressions for regular grids by using the expressions for 
unstructured grids and including the known areas, volumes, and cell face orientations. 
Other types of grid systems such as generalized curvilinear grids may also be viewed as 
specializations of the unstructured grid approach. 
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10.3 A Multigrid Scheme 


Implicit methods as currently practiced require dimensional splitting (approximate 
factorization) to be practical. For this reason, the unstructured meshes of the previous 
section, are not conducive to implicit methods. Because of CFL restrictions, explicit 
methods can be prohibitively expensive. An alternative of current interest is the use of 
multigrid methods. 

If any one feature characterizes multigrid schemes, it is the way they capture the gross 
features of a solution using a relatively coarse grid, using the usual fine grid to sharpen 
the solution. This approach is economical for two reasons. First, the coarse grid, having 
fewer unknowns, requires less arithmetic in each time step. Second, and just as important, 
those time steps can be bigger for a given CFL number. Looked at another way, if shocks 
can only move one grid point per step (a common restriction on explicit schemes), it helps 
to have fewer grid points to traverse. 

The second law provides an interesting insight into the stability of certain multigrid 
schemes. In particular, given a fine grid scheme which satisfies a semi-discrete cell entropy 
inequality, it is possible to construct multigrid schemes which also do. These multigrid 
schemes, termed unigrid schemes (ref. 39), were first introduced by Steve McCormick. 

In a unigrid scheme, each coarse grid cell consists of a collection of adjacent fine grid 
cells. Coarser grids contain more fine grid cells, but all grids are composed of a collection 
of cells from the finest grid. Consider the situation shown in figure 10.3. A single coarse 
grid cell is represented, which consists of four fine grid cells. 
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To advance in time on the fine grid is not difficult. Section 10.1 showed that schemes 
which satisfy a cell entropy inequality can be constructed. To advance in time on the 
coarse grid, consider the coarse grid cell as a single control volume. With this view, it is 
possible to advance in time without additional assumptions. The surface integral of the 
fluxes can be computed as on the fine grid using the fact that q is piecewise constant along 
the edges. In figure 10.3 the surface integral would be done in eight parts. 

For the moment, define qjK as the average value over the coarse grid cell. The semi- 
discrete conservation laws for cell JK can be written 

(qjK),t + ”(/j+|fc+i + /j+f * “ 

+ 2^(/jfc+f + /j+ifc+f ” fjk-\ ~ fj+ik-i) = 0 (10.3.1) 

where Ax and Ay refer to the fine grid spacing. The problem which remains is to apportion 
this change among the four cells which make up the larger control volume. The natural 
choice is to apportion the change according to the volume of each cell. This means 

(qjic),t = = (?jfc+i)jt = (fy+ifc)’* = (9j+ifc+i)»t (10.3.2) 


This choice may or may not satisfy the second law for the larger control volume. The 
second law may be written as 


(SjK),tvjK + Ay{F j+ i k+1 + F j+ | fc - -Fj-ifc+i “ F j-\k ) 

+ Az(F * + | + F j+lk+ 1 - F jk _i - F j+lk _ |) > 0 

(10.3.3) 

The first term, ( Sjx),tVjK represents an integration, carried out as follows: 

(SjK),tVJK =(Sjk),tVjk + (Sjk+l),tVjk+l 

+(*Sj+ifc),tVj+ifc + (S'j+i*+i),tt>j+ifc+i (10.3.4) 

where vjk is the volume of the coarse grid cell and Vjj,, Vjfc+i , Vj+i*, and are the 

volumes of the four fine grid cells. In this case all four are the same, just Ax Ay. The flux 
terms of (10.3.4) are simply the net export of entropy by convection, the area weighted 
sum of the entropy flux around the (in this case) eight segments of the control volume 
boundary. 

If (10.3.3) is not satisfied, the fault cannot lie with the flux terms. This is so because 
the second law can be satisfied in each cell separately, using the same fluxes. The only 
remaining suspect is ( Sjk ),t which reflects the way that increases in q are distributed. 
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The term ( Sjk )>t can be increased by mixing a fraction of the total fluid in each cell. 
The mixed fraction will have an increased entropy. If the mixed fraction reaches unity, the 
control volume entropy reaches a maximum. Since entropy production depends linearly 
on this fraction, it is possible to precisely compute the fraction required. A small amount 
of additional dissipation occurs when the mixed fluid is returned to the cells, and a cell 
average is taken to determine the final piecewise constant values. 

This idea can be quantified using equations such as 

qjk = « qjK + (l - a)qjk (10.3.5) 

where is obtained when qjk is updated using (10.3.2), and qjK is the average state 
in the coarse cell. It is also the maximum entropy state for the coarse grid cell. It is 
easily verified that (10.3.5), if applied to each fine grid cell with the same value of a, is 
conservative. 

The effect of a mixing operation such as this is to increase the entropy within the 
coarse grid cell. In each fine grid cell 

Sjk > aS(qjK) + (1 - «)Sjfc (10.3.6) 

which can also be written 

Sjk ~ Sjk > a(S(qjK) - Sjk) (10.3.7) 

Summing these contributions from each of the control volumes leads to 

vjk(Sjk - Sjk ) > « Vjk(S(qjK) ~ Sjk) 

+Vjk+i(S(qjK) ~ Sjk+i) 

+Vj+ik(S(qjK) ~ Sj+ik) 

+Vj+ik+i(S(qjK) ~ Sj+ik+i) (10.3.8) 

where Sjk is the volume-weighted average of the entropy in the four fine grid cells, and 
Sjk is Sjk when (qjK),t is distributed according to (10.3.2). The value of Sjk that is 
required can be deduced from the second law. This leaves a as the only unknown. Treating 
(10.3.8) as an equality, one can solve for a. This value can then be used in (10.3.5) and 
the resulting distribution will satisfy a (coarse) cell entropy inequality. Following this 
procedure ensures that, coarse grid or fine, the domain entropy will not decrease. Thus, 
the stability conjecture of Chapter 2 applies. 

As steady state is approached (steady on the scale of the grid at least) the fluxes 
more nearly balance, and a becomes very small, vanishing at convergence. Away from 
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convergence however, shocks can propagate much more quickly on the coarse grid than on 
the fine grid. 

There is no particular reason why the coarse grid cells are limited to four fine grid 
cells. The only requirement is that the coarse grid must use the same volume and flux 
integrations that the fine grid uses. This ensures that convergence on the finest grid implies 
convergence on all grids. 


10.4 Algorithm Improvements 

None of the schemes shown in this work exhibit more than second-order accuracy. 
There seems to be a need for more accurate schemes, particularly when simulating tur- 
bulence without algebraic modeling. This translates into schemes which have less entropy 
production, for the case of purely convective problems. There are severed ways to modify 
the schemes shown so far, to reduce the cell entropy production rate. 

Let’s briefly review the second-order schemes shown in Chapters 5, 6, and 7 These 
methods consist of three parts. First, a “target” value of Qj+i is chosen. Second, the 

entropy production rate in each cell (or half cell) is computed. Finally the value of <lj+i is 

modified as necessary to satisfy a cell entropy inequality. To achieve higher accuracy, one 
or more of these steps needs to be improved. 

Perhaps the simplest thing to change is the target value. Up to now, it has been 
tacitly assumed that the “best” value for 1 (for scalar equations in one dimension) is 
the average value, |(ttj + Uj+i). This is the best from a Taylor series viewpoint, but there 
are other viewpoints. 

Recall that 

(P.)j =(>. + ); + (*>.-); (10.4.1) 

One constraint which is fairly easy to satisfy (in principle) is 

(P. + )i=-(Oj+i (10-4.2) 

Both sides of the equation depend only on 1 . If Uj + x is chosen to satisfy (10.4.2) then 
the global entropy production rate will vanish, as for Tadmor’s scheme. As discussed in 
Section 9.3 and 9.4, this property guarantees stability but can lead to unphysical solutions. 
Deviations from the target value will probably be necessary to satisfy a local entropy 
inequality. On the other hand, as the mesh becomes finer and finer, these deviations 
become smaller and smaller (for smooth solutions). This means that unphysical global 
entropy production vanishes in the limit as the mesh is refined, a desirable property. 

For illustration, consider the linear wave equation again. It happens that the value 
of Uj + i which satisfies (10.4.2) is the average value. If the target value is used without 
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modification, the resulting scheme is algebraically equivalent to central differencing with 
a Crank-Nicolson time advance. Such a scheme is known to be unconditionally stable but 
far from TVD, generating severe oscillations. 

When the values of u J+ i are modified to satisfy a cell entropy inequality as well, 
the resulting scheme is the second-order accurate one discussed in Section 5.7. Typical 
solutions are shown in figure 5.8. Thus, (10.4.2) has already been tried for the wave 
equation. It has not led to any schemes with better than second-order accuracy. While 
simple and elegant, (10.4.2) is not enough by itself. 

Another way to pick target values is to include more points in the stencil. Perhaps 
the ultimate in schemes of this type is the subcell resolution method of arbitrary order 
proposed by Harten (ref. 40), and Harten, Enquist, and Osher mref41. Although such 
schemes are essentially nonoscillatory (ENO), they are known to contain small violations 
of local entropy conditions. The gradient following techniques discussed in Section 5 show 
how to modify u J+ i so as to correct these small defects. Unfortunately, the ENO approach 
is restricted to one dimension. 

Regardless of which target value is selected, the method outlined in Section 5.3 can be 
used to modify it in such a way that a cell entropy inequality is satisfied. Since this involves 
finding a root of a nonlinear equation, iteration is generally called for. Unfortunately, after 
all the work of solving for <j> to several decimal places, the resulting dissipation is usually 
about twice what is required (see fig. 7.2). This is because the increase in (P+ )j is about 
equal to the increase in (P“)j+ i (at <f> = 1 the slopes are equal if entropy variables are 
used). To first-order, the increase in (P“)j will be equal to the increase in (P+ )j. In the 
worst case, however, (P“)j doesn’t change, and (P+)j must allow for this. The double 
correction occurs because of this disparity between the worst case and the usual case. 

This suggests under-relaxing the first nonlinear iteration by a factor of two and re- 
calculating (P*)j. After this iteration (P»)j will be very close to zero in most the cells in 
which it had been negative. It is recommended that this underrelaxation be discontinued 
after one iteration. 

A third way to reduce entropy production is to extend the method suggested in Section 
5.3 in the following way. After u J+ 1 has been modified to make (P t )j > 0 in each cell, (P t )j 
can be reevaluated. This value can be used to give each half cell an “entropy budget” with 
which to move <f>j+i back toward the target value, thereby reducing the entropy production 
in cells j and j + 1. The variations on this third approach are nearly infinite, but they are 
bound to be costly. 


10.5 Grid Refinement 


One of the things which makes CFD an art instead of a science is the construction of 
a suitable computational mesh. It is often very difficult to construct a regular mesh about 
a complicated geometry. Small changes in geometry entail large changes in such grids. As 
the complexity increases, control over placement of grid points decreases dramatically. By 
contrast, construction of an unstructured mesh is not at all difficult, especially if mesh r 

clustering is not an issue. 

A promising approach is to produce a solution dependent clustering by moving the * 

existing grid points (ref. 42), or by adding new ones (ref. 43). With these two techniques, 
a computational mesh can be modified to improve resolution of steady or transient flow 
features of interest. Some researchers have gone so fax as to generate a new mesh every 
time step (ref. 44). 

Two problems have always hampered this approach. The first is the selection of a 
suitable sensor for detecting underresolved and over resolved areas. The second problem 
is the tendency to place all the points in the shock, which in the inviscid case is a true 
discontinuity and therefore unresolvable. 

The second law provides an elegant sensor, at least for inviscid problems. Such prob- 
lems produce no entropy except at shocks. This can be shown analytically. Numerically 
though, some (hopefully positive) rate of entropy production exists in each cell. Ignoring 
shocks for the moment, this cell entropy production rate (P t )j is a measure of resolution. 

The higher (P t )j is, the more poorly resolved the solution is locally. 

The problem of course is shocks. Shocks have enormous entropy production rates 
(both physically and numerically) and tend to invalidate the sensor. Standard techniques 
include imposing a minimum cell size (for mesh movement techniques) or simply refin- 
ing all cells (including those with shocks) which have more than some threshold entropy 
production rate (for mesh refinement techniques). 

10.6 Navier-Stokes Equations 

In many flows of interest, it is important to account for the effects of viscosity and 
heat transfer. The inclusion of these effects is what separates the Navier-Stokes equations 
from the Euler equations. This section will explore the effects of viscosity and heat transfer 
on the second law. 

The viscous terms are stresses, forces that act on the boundaries of each control * 

volume. These can affect the momentum and kinetic energy in the volume, which in turn 
influence the rate of entropy accumulation, S, t - They do not, however, directly affect the 
entropy flux F because they do not change the state of the fluid at the cell boundary. 


132 



The major effect of viscous stresses acting on a cell boundary is to reduce velocity 
differences between adjacent cells. Since momentum is conserved, this implies a reduction 
in the total kinetic energy of the two adjacent cells. Since energy is conserved the lost 
kinetic energy must appear as heat. The process of converting kinetic energy into heat 
produces entropy. 

The heat transfer terms appear in the energy equation and thereby influence the rate 
of entropy accumulation S,t . They are part of the definition of entropy flux. In particular, 
the term — ^ is added to the entropy flux F. The heat flux Q is usually determined by 
Fourier’s law; it is proportional to (but of opposite sign to) the temperature gradient. 

It is well known that entropy is always produced by heat transfer over a finite tem- 
perature difference. Given a heat flux Q between two cells of different temperatures, the 
entropy flux leaving the high temperature cell is ^ while the entropy flux entering the 

low temperature cell is ^ . Both temperatures are measured in degrees Kelvin or some 
other absolute temperature scale. Since the heat flux is generally from the hot cell to the 
cold cell, the entropy of the cold cell increases by more than the entropy of the hot cell 
decreases. This difference is the entropy production due to heat transfer. 

The inclusion of viscosity and thermal conductivity has interesting consequences. For 
example, since some entropy is produced “legitimately,” it becomes possible to scale back 
the “artificial” dissipation because this was added only for stability reasons. In some cases 
with very fine grids or very low Reynolds numbers, it may be possible to do away with it 
altogether. An extreme example is the heat equation for which no artificial dissipation is 
required. 

Secondly, including these terms may prove to be a far easier way to introduce entropy 
production than the way shown in Chapter 7. The Reynolds number and Prandtl number 
can be used as locally adjustable parameters which are set in some way which satisfies the 
cell entropy inequality. Similar ideas have been proposed by Dulikravitch (ref. 45). 

10.7 Chemistry 

Recently there has been some interest in computing very high speed flows ( M > 10) 1 
At the temperatures involved, oxygen and nitrogen react in ways that would be extremely 
unlikely at room temperature. This necessitates a different sort of entropy accounting. 

The first major change is that the convective entropy flux (as opposed to the heat 
transfer term) must be computed separately for each species being tracked. That is, the 

1 It is debatable whether such flows can be described by a continuum approximation at all since the Knudsen 
number is too big. The mean free path can exceed the boundary layer thickness in some cases. This issue 
is ignored here. 
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entropy carried by a molecule of oxygen is different than that of a molecule of nitrogen. It 
is no longer enough to know the mass transferred, the mass flux must be known for each 
species individually. 

The second major effect is that reactions occur within the control volume at a rate 
that depends on the concentrations of the reactants in the control volume. These reactions 
produce entropy. This must be reflected in the integration of the entropy within each cell 
volume. Each species must be accounted for separately as in the flux computation. 

10.8 Stability and Unsteady Flow 

The stability analysis presented here is more suitable for unsteady flows than is the 
more traditional von Neumann analysis. It is also more suitable for the nonlinear dissipa- 
tion techniques currently practiced. 

Traditional von Neumann analysis is fairly simple and therefore quite popular. It 
consists of supposing that the solution consists of a sine wave, and then applying the 
numerical scheme to this solution. If it can be shown that the solution will not increase 
in amplitude, regardless of frequency, then we can say with confidence that the scheme is 
stable in sense of Lax. The L 2 norm in wave space is bounded for all time. 

There are several problems with this. First, it is generally difficult to apply von 
Neumann analysis to nonlinear numerical schemes. By their nature, nonlinear schemes do 
not maintain the frequency of a sine wave. Even if a single frequency is used as an initial 
condition, the solution after one step contains many frequencies. The use of nonlinear 
dissipation is important. For example, it is impossible to construct a TVD scheme (even 
for one-dimensional wave equations, where this property is appropriate) which is more 
than first-order accurate, without using nonlinear dissipation (ref. 29). 

In many physical flows, disturbances of certain frequencies actually grow. This is 
resolved ultimately by transfer of energy into frequencies which don’t grow. Von Neumann 
stability analysis, by requiring stability in each mode separately, does not allow accurate 
simulation of such flows. 

By contrast, the satisfaction of a cell entropy inequality is no harder for nonlinear 
schemes than for linear schemes. The transfer of energy between frequencies does not 
affect the analysis. In addition, there are no physically allowable flows which violate the 
second law. Thus, there is no danger of excluding an interesting physical flow on the 
grounds that it is unstable. 


10.9 Summary 

There is considerable distance between the proposals of this chapter and the reality 
of a working simulation. (This is why this dissertation contains no results in higher di- 
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mensions or for more complicated equation sets.) Nevertheless, the possibilities appear 
bright; much more can be done than was previously possible. Provided that the conjecture 
of Chapter 2 holds, Chapter 7 demonstrates a scheme with true nonlinear stability for 
the Euler equations in one dimension. It appears that such schemes are extendable to 
multi dimensions without dimensional splitting. Extension to the Navier-Stokes equations 
also seems possible. 

The basic principle, a cell entropy inequality, appears sound for any problem in con- 
tinuum mechanics. By contrast, standard upwind schemes are limited to equations which 
have real eigenvalues and linearly independent eigenvectors. For practical purposes the 
eigensystem must be expressed in closed form. This may prove to be a serious limitation. 

Many schemes in use today rely on a geometric or intuitive definition of smoothness. 
These are subject to numerous pitfalls, the most serious of which is the lack of any con- 
nection to the physics. Many of them are also limited to one-dimensional representations. 
The cell entropy inequality approach does not share these pitfalls. 

Where the upwind methods shine is in one dimension. Here the de facto similarity 
between schemes which satisfy a cell entropy inequality and upwind schemes already in 
use is striking. From this point of view, the second law is merely a way of generalizing 
the schemes already in use. This, and the fact that there are no obvious roadblocks to 
the really interesting problems, leads to the conclusion that the ideas presented here will 
provide fertile ground for future research. 
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Chapter 11. SUMMARY 


In the author’s view, CFD is presently an art with elements of science. There are 
several steps in which “a miracle occurs.” A particular strategy is advocated without good 
reason, simply because it seems to work and nothing more satisfying is available. This work 
has addressed one such area, that of artificial dissipation. Artificial dissipation is analyzed 
from the viewpoint of numerical satisfaction of the second law of thermodynamics. 

A strong conjecture of nonlinear stability is offered for schemes which satisfy the 
second law over the entire domain (the cell-by-cell entropy inequality is not required for 
this). This conjecture is currently limited to periodic boundary conditions. 

Subsequent investigation provided indications that a cell-by-cell application of the 
second law provides more satisfactory solutions. Part of the difficulty in providing a proof 
lies in defining what is meant by “satisfactory.” For example, consider solution of the scalar 
wave equation, using three-point central differencing and Crank-Nicolson time advance. 
Such a scheme produces no net entropy over the domain, the physically correct result. It 
does, however, have a negative entropy production rate in some cells and a positive rate 
in others. The solution is bounded as required by the proof (and by linear theory), but 
it contains n um erous oscillations which are not physically correct. On the other hand, if 
this scheme is modified slightly to require a positive entropy production rate in each cell, 
the oscillations do not appear. Similar observations were made for a nonlinear, entropy 
conserving scheme in Chapter 9. There it .was noted that the wiggles were correlated in 
time, space, and magnitude with the negative entropy production rate in certain cells. 

The schemes suggested by this approach were formally shown to have the property of 
diminishing total variation as defined by Harten when applied to the linear wave equation. 
In the case of the nonlinear wave equation, the TVD property could be achieved with 
the assumption of locally constant wave speed. For the Euler equations in one dimension, 
there are still some similarities, but one major difference appears: while it isn’t possible to 
construct a scheme which is TVD for this problem, it is possible to satisfy a cell entropy 
inequality. 

Satisfaction of a cell entropy inequality for the one-dimensional Euler equations was fa- 
cilitated by the use of a surprising identity (7.3.1). This was discovered using MACSYMA, 
not derived with great insight. The author is not aware of any previous mention of this 
identity in the literature. It was pointed out by Roe that this identity also applies to 
the shallow water equations. Such an identity provides an unambiguous scaling of the 
eigenvectors in such problems. This is helpful in theoretical work but not too important 
in practical schemes. 
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Various extensions appear possible including multidimensions, viscosity, heat transfer, 
and chemistry. This flexibility is possible because of the widespread applicability of the 
second law. Since dimensional splitting is not required, unstructured grids become feasible. 
These may be refined using the entropy production rate itself as an indicator of local 
resolution. 

In addition to artificial dissipation and grid refinement, topics covered in this work, 
there are several more “black boxes” which are used in CFD. These are turbulence modeling 
(or turbulence itself!), boundary conditions, and nonlinear accuracy estimates. Until these 
are also understood, Computational Fluid Dynamics will remain on a house of cards, 
depending for validation on the very experiments it is supposed to render unnecessary. 
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Appendix A. DERIVATION OF SUFFICIENT CONDITIONS 

FOR TVD SCHEMES 


The objective here is to show that the total variation of the solution will not increase 
if the sufficient conditions are met. 

The derivation begins with the assumption that the independent variable is a scalar 
quantity u and that it can be advanced in time on a discrete mesh using 

«;+■ = u" - A« i .j + ; (A.1.1) 

where C,_ 1 and D. , 1 can be functions of u at the grid points. 

The objective is to derive conditions such that 

TV(u n+1 ) < TV(u n ) (A. 1.2) 

where TV is the so-called total variation of u. Total variation is defined as a sum so that 
(A. 1.2) can be written as 


E i«"+? - “? +, i * E K+i - «*i 


(A. 1.3) 


j j 

The sum is carried out over all values of j , that is, over all mesh intervals. 


The next step is to eliminate those variables which are at time level n + 1. Using 
(A.1.1) it is possible to determine by a simple index shift. 


u 7+i = “”+i ~ C j+ iAu j+ i + D j+ 1 Au i+ | 
Subtracting (A.1.1) gives an expression for 


(A. 1.4) 


“i+i - « n+1 - ( ' ,n 


J =K+I - ) + Cj-iAuj.i 

~ ( G i+\ + D j+i) Au j+± + D i+\ Au i+\ 


Collecting terms gives 

*?+ 1 - u l +1 = c i-\ Au i-\ + ( 2 - c i+\ ~ D i+\) Au i+\ + D i+i Au i+i 


(A. 1.5) 


(A. 1.6) 


At this point one takes the absolute value of both sides and applies the sum inequality for 
absolute values 


K+i 1 - “7 +1 l < K_1 Au^jl + 1(1 - C iH - D i+i ) A« i+j | + |B i+ jA» y+j | (A. 1.7) 
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Summing both sides, 


i i 

+ i(i - c i+\ - D j+±) Au j+± i 

+|D, + ,A« j+ ,|) (A. 1.8) 

Splitting the products and rearranging terms gives 

TV(«*+ , )<X)(|C < _ 1 ||Aii i _j|) 

i 

i 

+ £(1^+111^+10 (A.1.9) 

3 

At this point the key observation is made. Notice that, for example 

£(l C i-ill A “i-*l) = £(|C J+1 ||A« i+ ,|) (A. 1.10) 

i i 

except for boundary conditions. The boundary conditions are ignored in this derivation 
so the indices j — \ and j + | can both be shifted to j + 1. This means that (A. 1.9) may 
be written as 

TV(u"«) <£(|C i+1 ||Auj +1 |) 

i 

+ E(l(l-O i+ 1 -JJ, + 1 )l|A« J+ j|) 

i 

+ SO- D i+ill A «i+il) (A. 1.10) 

j 

which allows a massive collection of terms 

TV(«*+ , )<E(l a J+il + K 1 - Cf l+l -^+j)l + l B l+il)|A«>+jl (A. 1.10) 

j 

The only thing which prevents further collapse of the right-hand side is the absolute value 
operator. If the three terms are individually positive the absolute value operator is redun- 
dant. That is if 


c Hi > 0 

J 2 0 

C iH + D Hi < 1 (A. 1.11) 


139 



for all j, then the term in parentheses evaluates to unity and 

TV(u n+1 ) < TV(u n ) (A.1.12) 

Conditions (A. 1.11), first derived by Harten, are known as the sufficient conditions 
for TVD schemes. 
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